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Abstract 

The  Air  Force  Institute  of  Technology  (AFIT),  working  to  meet  requirements  set 
by  the  Air  Force  Research  Laboratory’s  Next  Generation  Engine  (NGE)  initiative,  is 
developing  upper  stage  rocket  models.  The  current  path  of  investigation  focuses  on 
combining  a  dual  expander  cycle  with  an  aerospike  nozzle,  or  the  Dual  Expander 
Aerospike  Nozzle  (DEAN)  using  methane  fuel.  The  methane  DEAN  (MDEAN)  design 
process  will  rely  heavily  on  AFIT’s  previous  work,  which  focused  on  the  development  of 
tools  for  and  the  optimization  of  a  liquid  hydrogen/liquid  oxygen  DEAN  engine.  The 
work  outlined  in  this  paper  expands  the  existing  research  by  substituting  liquid  methane 
for  liquid  hydrogen.  The  targets  derived  from  the  NGE  program  include  a  vacuum 
specific  impulse  of  383  seconds,  25,000  lbf  of  thrust,  and  a  thrust  to  weight  ratio  of  108. 
Additionally,  although  no  quantifiable  requirements  are  provided,  reliability  and 
reusability  are  identified  as  design  priorities.  NASA's  Numerical  Propulsion  System 
Simulation  (NPSS™)  was  used  in  conjunction  with  Phoenix  Integration's  ModelCenter™ 
to  optimize  over  several  parameters  to  include  O/F  ratio,  thrust,  and  engine  geometry. 
After  thousands  of  iterations  over  the  design  space,  the  selected  MDEAN  engine  concept 
has  349  s  of  Isp  and  a  thrust  to  weight  ratio  of  120.  The  MDEAN  was  compared  to  liquid 
hydrogen  technology,  existing  methane  technology,  and  the  NGE  goals. 
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METHANE  DUAL  EXPANDER  AEROSPIKE  NOZZLE  ROCKET  ENGINE 


I.  Introduction 

The  US  Air  Force  is  increasingly  reliant  on  capabilities  provided  through 
platforms  located  in  the  space  domain.  Mission  areas  currently  include  functions  vital  to 
the  joint  force  such  as  intelligence,  surveillance  and  reconnaissance,  command  and 
control,  navigation,  weather,  counterspace,  and  communications  [1].  Spacelift,  the 
mission  area  by  which  the  AF  gains  access  to  the  space  domain,  relies  on  chemical 
rocketry.  Access  to  the  space  domain  is  predicated  on  achieving  a  change  in  velocity  of 
approximately  7.8  km/s  for  low  Earth  orbit  and  an  additional  change  in  velocity  specific 
to  the  desired  location  of  the  asset.  Chemical  rocketry  is  currently  the  only  viable 
technique  for  attaining  the  required  launch  velocity  change  and  therefore  warrants 
continued  study  in  search  of  increased  efficiency  and  reliability. 

Fiscal  constraints  necessitate  researching  reliable,  cost  effective  technologies  that 
can  meet  national  spacelift  requirements.  Though  significant  (on  the  order  of  $5,000  per 
lb)  [2],  the  cost  of  space  launch  is  a  relatively  small  fraction  of  overall  program 
acquisition  cost.  Therefore,  the  chief  spacelift  requirement  is  reliability.  Following  this 
logic,  overall  program  costs  must  decrease  given  financial  constraints,  but  the  effort  to 
save  program  funds  through  research  and  development  of  the  space  launch  mission  area 
must  not  decrease  the  likelihood  of  successful  launch. 

Requirements 

A  combined  government/industry  team  developed  targets  for  launch  capability 
performance  under  the  Next  Generation  Engine  (NGE)  program  with  the  overarching 
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goal  to  reduce  the  cost  of  space  launch.  The  program  strives  to  advance  the  perfonnance 
of  existing  technologies  to  better  support  launch  capability  and  thereby  decrease  the 
overall  cost  of  delivering  assets  to  space.  The  goals  of  the  NGE  program  important  to  the 
current  research  are  the  performance  parameters  pertaining  to  upper  stage  engines.  To 
represent  a  significant  improvement  over  the  state  of  the  art,  the  NGE  has  determined  an 
upper  staged  engine  will  need  a  specific  impulse,  Isp,  of  383  seconds,  25,000  pounds  (lbf) 
of  thrust,  and  a  thrust  to  weight  ratio  of  108.  Additionally,  reliability  and  reusability  are 
cited  as  requirements  without  specific  quantifiable  measures.  The  current  research 
employs  a  Dual  Expander  Aerospike  Nozzle  (DEAN)  concept  to  attempt  to  meet  these 
requirements.  It  is  important  to  note,  while  AFIT  has  looked  to  the  NGE  program  as  a 
source  of  research  goals,  AFIT  is  not  associated  with  the  other  participants  of  the  fonnal 
NGE  program  (financially  or  otherwise). 

The  DEAN  Model 

Past  work  has  occurred  at  AFIT  related  to  this  subject  matter.  Thus  far,  a  model 
was  developed  (see  Figure  1)  to  characterize  the  performance  of  a  hydrogen/oxygen 
DEAN  upper  stage  rocket.  The  model  was  then  upgraded  and  incorporated  into 
ModelCenter™  (a  commercial  modeling  tool)  in  an  effort  to  automate  parametric  study 
of  design  variables.  Another  upgrade  introduced  a  more  realistic  isentropic  contour  to  the 
geometry  and  further  refined  the  model.  Current  research  seeks  to  leverage  the 
exceptional  work  of  the  past  AFIT  team  to  create  a  more  efficient  hydro gen/oxy gen 
modeling  tool  and  subsequently  convert  it  to  operate  with  methane  chemistry  and  fluid 
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properties.  The  previous  efforts  will  be  described  in  detail  in  the  background  chapter  of 
this  paper. 


Figure  1:  DEAN  Schematic.  From  Hall,  Inspired  by  Simmons. 

The  primary  distinguishing  characteristic  of  the  new  rocket  design  as  compared  to 
prior  AFIT  research  is  the  replacement  of  cryogenic  hydrogen  fuel  with  cryogenic 
methane.  For  an  ideally  expanded  rocket  nozzle,  Isp  varies  directly  with  the  rocket 
characteristic  exhaust  velocity.  To  maximize  the  exhaust  velocity,  the  molecular  weight 
of  the  exhaust  products  can  be  minimized  through  propellant  selection,  and  because  of  its 
low  molecular  weight,  hydrogen  is  an  attractive  propellant.  Despite  the  performance 
advantage  of  hydrogen  over  methane,  launch  operations  are  complicated  through  the  use 
of  cryogenic  hydrogen.  Additionally,  the  foremost  source  (approximately  96%)  of 
hydrogen  is  the  processing  of  fossil  fuel.  Natural  gas  specifically  accounts  for  about  29% 
of  annual  hydrogen  gas  production  [3],  Advantages  of  cryogenic  methane  relative  to 
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other  options  are  examined  in  the  background  chapter  of  this  paper  and  perfonnance 
comparisons  are  made  in  the  results  chapter. 

Research  Objectives 

Current  research  will  support  the  NGE  upper  stage  methane  rocket  performance 
goals  through  the  following  objectives: 

1 .  Update  the  existing  hydrogen/liquid  oxygen  DEAN  Model  with  an  emphasis  on 
modularization  and  efficiency. 

2.  Update  the  existing  DEAN  Model  with  methane/liquid  oxygen  chemistry  and 
fluid  properties. 

3.  Perfonn  parametric  studies  on  the  updated  model. 

4.  Analyze  the  resulting  design  points  against  NGE  performance  goals  and  historical 
designs. 

Multiple  design  tools  will  be  employed  to  pursue  the  research  objectives. 
ModelCenter  by  Phoenix  Integration  will  control  the  input  and  results  from  multiple 
modeling  tools.  Additionally,  ModelCenter  has  the  capability  of  performing  statistically 
relevant  parametric  studies  of  the  DEAN  design  variables.  NASA’s  Numerical 
Propulsion  System  Simulation  (NPSS™)  will  calculate  parameters  for  the  various 
components  of  the  rocket  model.  Finally,  NASA’s  Chemical  Equilibrium  with 
Applications  (CEA™)  will  be  used  for  chemical  data  including  fluid  properties. 
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II.  Literature  Review 


Background 

The  fundamental  goal  of  this  research  is  the  development  of  a  model  that 
calculates  perfonnance  parameters  of  a  methane  rocket  employing  a  dual  expander 
engine  cycle  and  an  aerospike  nozzle.  It  will  be  shown  whether  performance  peculiar  to 
the  chosen  initial  design  options  is  appealing  for  continued  development,  provided  there 
is  an  advantage  in  performance  over  existing  concepts.  This  chapter  examines  rocket 
engine  theory,  previous  DEAN  research,  the  preliminary  design  decisions,  and  the  basis 
for  perfonnance  comparison. 

Rocket  Theory 

Ideal  Rocket  Equation 

The  basis  of  the  key  performance  parameter  (Isp)  is  derived  from  the  transfer  of 
momentum  from  the  rocket  to  the  rocket  exhaust  gases  [4].  To  derive  the  relationship 
between  the  task  to  be  performed  by  the  rocket  (a  change  in  velocity)  and  the  change  in 
mass  of  the  rocket  due  to  the  consumption  of  propellant,  external  forces  are  neglected  and 
the  total  momentum  of  the  system  is  assumed  to  be  static: 

HP 

LLr  system  _  _  ,  j  + 

dt 

This  suggests  that  any  momentum  that  is  removed  from  the  rocket  body  is  fully 
transferred  to  the  exhaust,  or  the  change  in  the  momentum  of  the  rocket  is  equal  to  the 
opposite  the  change  in  momentum  of  the  exhaust. 

(m  —  dm)dv  =  —dmve  (2) 
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The  variable  ve  is  the  exhaust  velocity  of  the  propellant  exiting  the  rocket,  and  it 
defines  the  primary  rocket  performance  parameter  through  the  following  relation. 


^spdo  (3) 

Simplifying  these  expressions  results  in  an  integral  where  mj  and  mf  are  the  initial 


and  final  mass  of  the  rocket,  respectfully. 


fAv  Cmf  dm 

dv  =  —vP  - 

Jo  Jm,  m 


(4) 


This  reduces  to 


Av  =  —veln  | 


e) 

The  final  simplification  is  the  substitution  relating  exit  exhaust  velocity  to  specific 


impulse: 


(mf\ 

Av=  -Ispg0ln^—J  (6) 

Now,  given  a  specific  mission  (Av),  a  rocket  (Isp)  can  be  chosen  that  will  provide 

the  needed  capability  with  an  estimate  of  the  necessary  fuel  (ny  -  mf).  The  challenge 
undertaken  in  this  study  is  to  calculate  a  reasonable  estimate  of  Isp  for  a  specific  rocket 
concept  over  a  range  of  design  variables. 

Possible  Mission  Scenario  and  Requirements 

Possible  mission  requirements  are  described  by  the  NGE  solicitation  for  an  upper 
stage  hydrogen  rocket  for  the  future  replacement  of  the  RL-10  [5].  The  most  basic  rocket 
mission  requirement  is  the  translation  of  the  mission  into  a  quantifiable  capability  that 
leads  to  engineering  decisions.  In  the  case  of  an  upper  stage  engine,  a  suitable  starting 
point  is  the  required  Av.  For  this  effort,  a  common  upper  stage  maneuver,  a  transfer  from 
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a  parking  orbit  (v;  =  7.79  km/s)  to  a  geosynchronous  orbit  (vf  =  3.08  km/s),  is  taken  as  a 
baseline  for  a  relevant  Av.  As  rocket  propulsion  has  relatively  high  thrust,  an  impulsive 
transfer  between  the  two  orbits  of  interest  with  an  inclination  change  of  28°  at  the  apogee 
of  the  transfer  orbit  is  used  to  calculate  the  Av.  The  inclination  change  is  included  based 
on  launching  out  of  Cape  Canaveral,  FL.  The  first  maneuver  is  calculated  as  a  Hohmann 
transfer  of  2.46  km/s.  To  calculate  the  scalar  change  in  the  velocity  of  the  two 
maneuvers,  the  law  of  cosines  is  employed  for  the  second  maneuver. 

Avb  —  J (vf  +  vj  —  2vtVf  cos  6?)  (7) 

Avtotal  =  4.29(km/s)  =  Ava  +  V(l-592  +  3.082  -  2(1.59)(3.08)  cos(28))  (8) 

The  overall  mission  associated  with  the  upper  stage  must  achieve  a  change  in 
velocity  of  4.29  km/s.  The  change  in  velocity  will  be  used  with  the  ideal  rocket  equation 
to  determine  if  an  engine  with  a  given  Isp  is  suitable  for  the  upper  stage  representative 
mission.  A  challenge  with  the  conceptual  development  of  useful  rocket  engines  is  the 
nature  of  the  space  enterprise.  Because  of  the  expense  and  the  global  reach  of  space 
assets,  many  payloads  are  unique.  Each  payload  needs  mission  analysis  to  detennine  the 
proper  propulsion  solution,  but  the  DEAN  model  development  assumes  a  general 
mission.  Although  the  current  work  doesn’t  fully  support  a  rigorous  systems  engineering 
process,  the  goal  is  to  develop  a  robust  enough  engine  concept  that  the  design  will  satisfy 
several  different  space  lift  missions,  or  at  the  very  least,  find  an  interesting  niche  for 
which  it  is  particularly  suitable. 

The  NGE  solicitation  provides  additional  requirements  in  the  form  of  constraints 
for  an  upper  stage  hydrogen  engine.  Total  thrust,  F,  is  to  be  set  at  25,000-35,000  lbf  [5]. 
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Additionally,  the  thrust  to  weight  ratio,  T/W,  is  to  be  greater  than  or  equal  to  108,  and  Isp 
must  be  greater  than  465  seconds.  While  the  NGE  solicitation  specifically  identifies  a 
hydrogen/oxygen  cryogenic  upper  stage  engine,  it  is  still  of  interest  to  investigate  the 
possibility  of  a  methane  rocket.  NGE  seeks  methods  to  increase  capability  over  a 
baseline,  the  RL10.  The  current  work  derived  derives  the  MDEAN  requirements  from 
similar  increases  in  perfonnance  over  a  different  baseline,  the  RD-185,  while  keeping  the 
basic  functionality  of  providing  25,000  lbf  of  vacuum  thrust.  The  RD-185  has  a  vacuum 
Isp  of  378  and  a  thrust  to  weight  ratio  of  44  [6].  These  baseline  requirements  will  be  used 
as  input  variables  for  parametric  study,  constraints  for  the  design  (thrust),  or  calculated 
output  (Isp,  T/W).  The  NGE  requirements  are  the  criteria  used  to  evaluate  the 
perfonnance  of  the  MDEAN. 

Engine  Cycle 

The  characterization  of  the  method  of  propellant  transfer  into  the  combustion 
chamber  is  known  as  the  engine  cycle  (excluding  the  specifics  of  injection).  Several 
engine  cycle  options  exist  for  a  cryogenic  liquid  rockets.  The  driving  factors  in  the 
refinement  of  cycle  selection  are  the  range  of  pressurization  necessary  from  the 
propellant  tanks  to  the  combustion  chamber,  geometry,  complexity,  sizing,  and 
perfonnance.  Three  engine  cycles  are  more  popular  and  have  been  demonstrated  in 
mission  capable  designs.  Diagrams  of  the  three  most  popular  cycles  (gas  generator, 
staged-combustion,  and  expander)  can  be  found  in  Figure  2. 
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Figure  2:  Common  Rocket  Cycles. 


One  top  level  distinction  is  the  open  system  versus  the  closed  system.  The  open 
cycle  combusts  a  portion  of  the  propellant  flow  to  drive  a  turbine  which  in  turn  operates 
the  fuel  and  oxidizer  pumps.  The  products  from  the  turbine  are  exhausted  overboard,  and 
this  cycle  concept  is  known  as  the  gas-generator  cycle.  This  cycle  provides  a 
comparatively  simple  solution  to  propellant  pressurization  at  the  cost  of  perfonnance. 

For  a  given  chamber  pressure,  the  gas-generator  cycle  decreases  overall  engine  Isp  by  2% 
to  5%  [4].  The  driving  force  for  concept  refinement  in  the  early  stages  of  DEAN  was 
based  on  maximizing  Isp  at  a  specific  thrust  and  thrust  to  weight  ratio.  Because  the  goal 
of  DEAN  is  high  performance,  the  gas-generator  cycle  was  not  selected  for  modeling. 

Closed  cycles  have  greater  opportunity  to  reach  the  high  level  of  performance 
required  by  the  NGE  solicitation  [4].  The  two  most  common  closed  cycle  concepts  that 
are  examined  are  the  staged-combustion  and  the  expander  cycles. 
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The  expander  cycle  flows  propellants  through  a  heat  exchanger  which  cools  the 
combustion  chamber  and  nozzle.  The  coolant  is  traditionally  the  fuel  which  absorbs  heat 
and  drives  a  turbine  which  operates  the  fuel  and  oxidizer  pumps.  All  of  the  unreacted 
propellant  then  enters  the  combustion  chamber  and  is  expanded  through  the  nozzle.  The 
DEAN  design  concept  differs  from  the  traditional  expander  cycle  in  its  dual  nature.  The 
DEAN  flows  both  the  fuel  and  oxidizer  through  heat  exchangers  and  each  flow 
independently  powers  the  fuel  and  oxidizer  pumps.  Overall  system  weight  is  reduced 
when  compared  to  the  other  cycles  because  there  is  no  gas  generator  or  extra  combustion 
chamber  providing  energy  to  the  turbo  machinery.  Additionally,  the  organization  of  the 
regenerative  cooling  system  associated  with  the  cycle  has  a  synergistic  match  with  the 
geometry  inherent  to  an  aerospike  layout.  The  layout  of  the  engine  plumbing  is  therefore 
simplified  to  an  extent,  and  the  likelihood  of  fuel  and  oxidizer  interaction  outside  the 
combustion  chamber  is  reduced.  The  probability  of  this  known  failure  mode  is  therefore 
reduced,  theoretically  increasing  the  reliability  of  the  overall  rocket  system. 

The  staged  combustion  cycle,  as  seen  in  Figure  2,  has  a  similar  flow  to  the 
expander  cycle  in  that  the  propellants  travel  through  a  heat  exchanger  which  is  linked  to 
the  combustion  chamber  and  nozzle.  After  the  propellants  are  heated,  a  portion  flows 
into  a  pre-bumer  combustion  chamber  that  consumes  the  flow  of  fuel  and  some  of  the 
flow  of  oxidizer.  The  products  of  this  combustion  operate  a  turbine  to  power  the 
propellant  pumps  and  then  flow  into  the  primary  combustion  chamber.  The  primary 
advantage  of  the  staged  combustion  concept  is  the  high  achievable  chamber  pressure 
which  will  lead  to  a  high  performance  engine  [7].  The  primary  disadvantage  is  the 
volume,  complexity,  and  mass  of  the  extra  combustion  equipment.  Although  this  cycle 
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would  conceivably  work  for  the  DEAN,  the  chosen  mission  scenario  of  upper  stage  space 
launch  favors  a  smaller,  simpler  concept.  The  ideal  design  concept  will  be  small  enough 
to  fit  within  the  confines  defined  by  existing  hardware,  not  to  exceed  90  inches  in  length 
with  a  diameter  not  to  exceed  73  inches  [5]. 

Therefore,  having  eliminated  the  perfonnance  drop  of  the  gas  generator  and  the 
bulk  of  the  staged  combustion  cycles,  the  DEAN  concept  incorporates  an  expander  cycle. 
The  additional  design  feature  of  the  dual  expander  system  takes  advantage  of  the 
aerospike  geometry  and  will  ease  plumbing  design. 

Turbo  machinery 

Turbo  machinery  refers  to  the  turbines  and  pumps  that  make  up  the  expander 
cycle  of  the  DEAN.  This  system  of  equipment  is  highly  tuned  to  a  specific  rocket  design 
and  provides  propellants  at  the  proper  pressure,  temperature,  and  mixture  ratio  to  drive 
the  combustion  process  in  the  chamber  at  the  design  flow  rate.  The  primary  advantages 
of  using  turbopumps  in  a  rocket  engine  concept  is  their  capability  to  provide  high 
combustion  chamber  pressures  with  lower  propellant  storage  tank  pressures.  Both  of 
these  conditions  tend  to  decrease  the  overall  mass  of  the  engine  concept.  Lower  pressure 
tanks  can  have  thinner  walls  due  to  lower  hoop  stress.  A  higher  chamber  pressure  can 
equate  to  a  volumetrically  smaller,  potentially  lighter  chamber  with  improved 
perfonnance. 

Pressure  Budget 

To  begin  the  sizing  process  for  the  turbopumps,  the  pressure  drop  across  the 
system  must  be  estimated  [4].  The  estimation  process  outlined  by  Humble  et  al.  starts 
with  the  desired  chamber  pressure  and  steps  through  the  pressure  losses  through  the 
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system  in  reverse  order:  chamber  pressure,  injector,  turbine,  cooling  jacket,  pumps, 
propellant  tanks.  For  an  initial  estimate,  the  pressure  drop  across  the  injector  is  estimated 
at  20%  of  chamber  pressure  at  nonnal  operating  conditions  or  30%  of  chamber  pressure 
for  throttled  operation.  This  pressure  drop  is  due  to  one  of  the  injector’s  primary 
functions,  the  isolation  of  propellant  flow  perturbations  between  the  chamber  and  the 
feed  system  [4].  The  next  step  is  to  assume  a  pressure  ratio  across  the  turbine.  Humble 
et  al.  estimates  a  turbine  pressure  ratio  of  1.5  for  an  expander  cycle  [4],  Hall  estimated  a 
pressure  ratio  of  1.84  for  the  oxidizer  and  a  pressure  ratio  of  1.56  for  the  hydrogen  fuel 
[8].  Based  on  historical  data,  Humble  et  al.  recommend  a  cooling  jacket  pressure  drop  of 
15%  of  chamber  pressure  [4],  The  associated  dynamic  pressure  drop  can  be  estimated  by 
Bernoulli’s  equation: 


A p  -  —  pv2 

The  velocity  can  be  assumed  to  be  about  10  m/s  and  a  piping  diameter  can  be 
based  on  equation  10,  calculated  from  the  pipe  cross  sectional  area,  Apipe.  As  with 
equation  9,  p  is  density  and  v  is  velocity. 


(9) 


171  pv  Apipg  (10) 

By  assuming  a  reasonable  storage  tank  pressure,  the  required  pressure  increase 
due  to  pumping  can  be  calculated.  A  ratio  of  note,  the  pressure  ratio,  is  the  required 
pressure  entering  the  cooling  jacket  divided  by  the  storage  tank  pressure.  This  ratio  will 
be  used  in  an  engine  balance. 
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Engine  Balance 


The  pressure  rise  provided  by  the  pumps  corresponds  to  a  power.  The  energy 
imparted  to  the  flow  through  the  pumps  must  be  balanced  through  the  energy  imparted  to 
the  turbine  and  through  the  thennal  energy  captured  through  the  cooling  jacket.  The 
balance  starts  by  setting  the  required  pump  power  equal  to  the  power  of  the  turbine  [4] 
while  assuming  a  constant  specific  heat  over  the  range  of  encountered  temperatures: 

g0mHp 

Vreq  ~  ~  —  VT^-T^p^i 

Vp 

where 

Preq  =  Power  required  to  drive  the  pump 
g o  =  Acceleration  due  to  gravity  (m/s') 
m  =  Mass  flow  rate  (kg/s) 

Hp  =  Pump  head  rise  (m) 
r|p  =  Pump  efficiency 
r|x  =  Turbine  efficiency 
Cp  =  Constant  pressure  specific  heat  (J/kgK) 
y  =  Isentropic  parameter 
Tj  =  Turbine  inlet  temperature  (K) 

Ptrat  =  Turbine  pressure  ratio 

For  the  DEAN  concept,  each  turbine  must  be  balanced  against  the  corresponding 
pump(s).  The  turbine  inlet  temperature  must  be  detennined  based  on  the  energy  imparted 
to  the  propellants  through  the  walls  of  the  cooling  jacket.  The  solution  must  be  found 
through  iteration  until  the  powers  match  to  within  acceptable  tolerance. 
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Existing  Work 

The  existing  design  provided  by  Martin  includes  two  turbopump  assemblies 
which  were  developed  by  Arguello  [9]  and  Strain  [10]  specifically  for  the  DEAN. 

Because  the  development  of  turbomachinery  specifications  is  laborious,  this 
research  leaves  specific  component  selection  and/or  development  to  future  research. 
Based  on  model  development,  reasonable  pressure  ratios  will  be  assumed  for  the 
turbomachinery  and  specific  calculated  changes  in  propellant  state  will  be  checked  for 
viability.  Past  work  has  developed  calculated  efficiencies  for  liquid  oxygen  and  liquid 
hydrogen  turbomachinery.  The  change  over  from  hydrogen  to  methane  will  significantly 
change  the  design  of  the  turbopump.  The  fuel  with  higher  density,  methane,  will  require 
less  head  and  therefore  less  power  [7].  This  will  reduce  the  size  of  the  pump  assembly 
for  the  methane  fueled  rocket,  saving  mass. 

Combustion  Chamber 

The  combustion  (or  thrust)  chamber  is  the  heart  of  the  rocket  where  the  chemical 
energy  of  the  propellants  is  converted  into  thrust  (see  Figure  3).  The  propellants  are 
atomized  through  the  injector  so  the  resulting  droplets  are  quickly  vaporized.  Then  the 
combustion  of  the  gaseous  propellants  increases  the  temperature  and  the  flow  rate.  The 
combustion  products  flow  is  then  accelerated  through  a  nozzle  to  supersonic  velocities 
prior  to  ejection  [11].  The  process  heats  and  accelerates  the  propellants  and  can  generate 
an  extreme  environment  of  high  pressures  and  temperatures  which  the  combustion 
chamber  must  survive  for  the  duration  of  the  mission.  The  added  qualitative  goals  of 
reusability  and  reliability  increase  the  need  for  the  chamber  to  resist  deformation  or 
failure  during  operation  over  multiple  cycles.  Increased  resilience  can  be  attained 
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through  a  higher  margin  of  safety  for  load  calculations  which  will  result  in  thicker  walls 
or  exploration  into  stronger  materials. 


Figure  3:  Thrust  Chamber. 


Performance  Measures 

An  important  quantitative  measure  of  a  thrust  chamber  is  the  characteristic 
velocity,  c*.  It  will  be  shown  later  (equations  26  through  28)  that  c*  is  proportional  to  the 
primary  rocket  measure  of  performance,  Isp. 


P<A  Vc*^YRTc 


(12) 


sfc 

c  is  dependent  on  the  components  of  the  product  gas  at  the  exit  of  the  chamber. 
Once  a  propellant  has  been  chosen,  every  variable  will  be  set  (within  a  range)  in  the 
above  equation  except  for  the  temperature.  The  characteristic  velocity  then  mostly 
depends  on  the  temperature,  which  should  be  maximized  to  the  extent  feasible.  Too  high 
a  temperature  will  weaken  the  chamber  walls,  causing  failure  due  to  the  extreme  forces 
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•  •  (  # 
involved  with  the  high  pressures.  For  methane  and  oxygen,  an  estimated  c  to  start 

engineering  estimates  is  6020  ft/s  [7]. 

Another  important  measure  of  the  chamber-nozzle  pair  is  the  thrust  coefficient, 
Cf,  which  represents  the  gas  expansion  performance  through  the  nozzle  and  is  defined  as 
the  force  generated  with  the  expansion  through  the  nozzle  over  the  force  due  to  the 
pressure  over  the  area  of  the  throat  alone,  or: 

F 

c,  =  JFc  <13> 

where 

F  =  Thrust  (lbf) 

At  =  Throat  Area  (in") 

Pc  =  Chamber  Pressure  (psia) 

Combustion  chamber  sizing 

One  possible  goal  in  designing  the  layout  of  the  combustion  chamber  is  to  ensure 
complete  combustion  prior  to  acceleration  through  the  nozzle  portion  of  the  chamber. 
Any  unreacted  propellant  exiting  the  combustion  chamber  represents  a  loss  of  potential 
energy  and  a  decrease  in  the  conversion  of  chemical  energy  into  thrust.  This 
phenomenon  is  complicated  by  the  possibility  of  incomplete  combustion  providing  a 
lower  average  molecular  weight  of  the  combustion  products,  which  also  increases 
perfonnance.  The  characteristic  length  is  related  to  the  completeness  of  combustion  in  a 
chamber  and  is  defined  as: 

_Ye_ 

~  At  (14) 

where 


16 


L  =  Characteristic  Length  (in) 

•5 

Vc  =  Chamber  Volume  (in  ) 

At  =  Throat  Area  (in“) 

^  .  •  •  • 

L  is  the  length  a  chamber  (of  equal  volume)  would  have  if  it  was  a  straight  tube 

without  the  converging  nozzle  section.  An  early  estimate  for  a  methane  rocket 

characteristic  length  can  be  derived  from  the  ranges  of  values  found  in  Humble  et  al  [4]. 

A  cryogenic  hydro  gen/oxygen  rocket  can  range  from  0.76  to  1.02  meters.  An  RP-1 

/oxygen  rocket  can  range  from  1.02  to  1.27  meters.  Larger  L*  can  cause  system  level 

perfonnance  degradations  by  increasing  size  and  weight  of  the  chamber,  increasing  the 

surface  area  requiring  cooling,  and  increasing  frictional  losses  in  the  chamber  [12]. 

The  residence  time  is  a  measure  of  the  mean  length  of  time  a  molecule  spends  in 

the  reactor  [12], 


PVC 


m 


(15) 


L  relates  to  the  residence  time  through 


v  =  Yl  =  WVctR 


(16) 


At  At 

The  residence  time  depends  on  the  chemical  kinetics  of  the  reacting  species  and  is 
detennined  experimentally.  The  residence  time  can  vary  between  0.0001  and  0.040 
seconds  for  different  types  and  sizes  of  chamber  [7].  As  the  characteristic  length  of  a 
chamber  increases,  the  volume  of  the  chamber  will  increase,  the  velocity  of  the 
propellants  will  decrease  and  the  residence  time  will  increase.  There  is  an  opportunity 
for  exploration  of  the  thenno  chemistry  with  this  relationship  and  L*  is  therefore  an  input 
design  variable  for  the  MDEAN.  For  the  MDEAN  modeling  effort,  the  L  design 
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variable  exploration  will  begin  in  the  neighborhood  of  one  meter,  the  high  end  for  a 
hydrogen  rocket  and  the  low  end  of  an  RP-1  rocket. 

Nozzle 

The  performance  measures  most  important  to  the  NGE  program  are  specific 
impulse  (Isp),  thrust  (F),  and  the  thrust  to  weight  ratio  (T/W).  While  the  rocket  as  a  whole 
must  be  considered  to  calculate  these  values,  the  equations  used  to  calculate  them  are 
largely  centered  on  the  characteristics  and  states  defined  by  the  chamber  and  the  nozzle. 
The  Isp  was  defined  with  the  development  of  the  ideal  rocket  equation  as  being  directly 
proportional  to  the  velocity  of  the  combustion  products  at  the  exit  of  the  nozzle.  Another 
representation  is  as  a  ratio  of  the  thrust  to  the  propellant  mass  flow  rate,  or 


Isp 


mg0 


(17) 


This  equation  describes  Isp  as  a  measure  of  the  total  perfonnance  of  the  rocket.  It  is  the 
generated  kinetic  thrust  imparted  to  the  rocket  per  mass  flow  rate  of  the  propellant.  A 
useful  analogy  to  this  relation  would  be  the  fuel  economy  of  an  automobile,  or  miles  per 
gallon. 


The  nozzle  accelerates  and  ejects  combustion  products  to  impart  momentum  to 
the  rocket.  Thrust  is  generated  by  a  rocket  through  two  mechanisms,  momentum  thrust 
and  pressure  thrust  [4],  Total  momentum  of  the  system  is  conserved  at  zero,  so  as  the 
fluid  accelerates  and  exits  the  rocket,  the  rocket  builds  momentum  equally  but  in  the 
opposite  direction  creating  momentum  thrust.  Pressure  thrust  is  created  as  the  pressure 
exiting  the  rocket,  pe,  can  be  unequal  to  the  ambient  pressure,  pa.  With  one-dimensional, 
steady  flow  through  the  nozzle,  there  will  only  be  a  pressure  differential  in  the  direction 
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of  propellant  flow.  Pressure  differentials  in  other  directions  will  be  equal  and  opposite 
due  to  the  symmetry  of  the  rocket  and  therefore  cancel.  The  basic  thrust  equation  then 
becomes: 


F  =  momentum  thrust  +  pressure  thrust  (18) 

F  =  mve  +  (pe  -  pa)Ae  (19) 

As  the  gas  exits  the  combustion  chamber,  it  enters  into  the  nozzle  and  is 

accelerated  to  Mach  1 .  After  passing  through  the  throat,  the  gas  is  further  accelerated 

through  the  divergent  portion  of  the  nozzle  and  decreases  in  temperature  and  pressure  as 

thennal  energy  is  converted  to  kinetic  energy.  The  process  over  the  length  of  the  nozzle 

is  complicated  as  the  thermodynamic  properties  are  dependent  on  the  constantly  changing 

state  of  the  gas.  Several  different  assumptions  can  aid  in  the  calculation  of  the  state 

through  the  nozzle. 

The  first  possible  simplification  for  the  flow  through  the  nozzle  is  the  condition  of 
frozen  composition.  The  assumption  asserts  that  the  composition  of  the  flow  remains 
constant  through  the  expansion  of  the  nozzle.  There  are  no  chemical  reactions  or  phase 
changes.  To  restate,  the  composition  at  the  exit  of  the  nozzle  is  identical  to  the 
composition  for  the  chamber  condition.  This  assumption  is  chosen  for  simplicity  of 
calculation  and  because  the  resulting  performance  estimates  are  conservative.  Frozen 
flow  calculations  underestimate  performance  by  one  to  four  percent  [11].  The  frozen 
flow  assumption  will  be  primarily  utilized  through  the  calculations  of  gaseous  properties 
provided  by  CEA. 
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Equilibrium  flow  is  more  complex  and  assumes  the  products  constantly  shift 
during  the  expansion  portion  of  the  nozzle.  The  equilibrium  is  modeled  to  exist  between 
chemical  composition  and  phase.  This  more  complex  assumption  results  in  inflated 
calculations  for  Isp  or  c*  and  is  therefore  not  used  in  the  DEAN  [11]. 

The  final,  unused  assumption  attempts  to  model  the  equilibrium  flow  to  include 
the  reaction  rates  of  the  steps  between  the  multiple  equilibrium  states.  The  complexity  of 
the  calculation  combined  with  the  lack  of  sufficient  data  on  the  rates  of  the  assumed 
reaction  mechanisms  usually  prevent  this  type  of  calculation. 

Aerospike 

The  aerospike  nozzle  has  several  advantages  over  a  conical  nozzle  and  was 
therefore  incorporated  in  the  initial  design  considerations.  The  development  of  the 
calculations  to  determine  aerospike  performance  is  explored  in  chapter  3. 

Conical  approximation 

For  the  purpose  of  this  study,  the  contour  of  the  aerospike  nozzle  will  be  treated 
as  conical.  Previous  DEAN  work  demonstrated  that  the  conical  aerospike  developed 
reasonable  approximations  for  mass  calculations  while  significantly  decreasing 
processing  time.  Specifically,  conical  approximations  vary  by  less  than  10%  [8]  from 
calculations  based  on  the  more  complex  contours.  Therefore,  to  save  computational  time 
and  operate  with  streamlined,  easily  testable  and  modifiable  code,  Hall’s 
recommendation  of  using  a  conical  nozzle  contour  will  be  implemented  in  the  current 
work. 
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Truncation 


Previous  DEAN  work  (the  first  DEAN  model)  explored  the  impact  of  truncating 
the  aerospike  nozzle  downstream  of  the  throat.  The  conclusion  was  that  Isp  didn’t  vary 
significantly  with  nozzle  length  and  the  truncation  provided  savings  to  the  mass  of  the 
engine.  A  truncation  of  25%  resulted  in  a  reduction  of  8.5%  of  the  mass  of  the  nozzle 
[13].  Reduction  of  the  supersonic  expansion  section  of  the  nozzle  of  up  to  75%  only 
reduced  modeled  Isp  by  1 .6%.  From  Martin’s  work,  it  can  be  concluded  that  a  modest 
truncation  of  25-50%  should  result  in  a  lighter  engine  with  slightly  lower  perfonnance. 


Figure  4:  Conical  Aerospike  vs.  Truncated  Conical  Aerospike. 


The  truncation  of  the  optimized  3ld  generation  DEAN  model  further  demonstrated 
that  within  a  regime  of  up  to  50%  truncation  (see  Figure  4),  performance  is  not 
significantly  reduced  [8].  Hall’s  model  estimated  performance  loss  of  approximately 
0.83%  for  49%  nozzle  truncation.  Truncation  beyond  this  length  impacted  the  expander 
cycle  portion  of  the  cooling  jacket  that  is  located  on  the  nozzle.  Significant  decreases  in 
performance  occur  beyond  80%  truncation.  Based  on  the  changes  in  performance  and 
mass,  Hall  concludes  any  change  in  the  truncation  of  the  nozzle  must  be  considered  in  the 
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context  of  a  full  stage  design  [8],  A  shorter,  lighter,  lower  perfonning  engine  might  be 
the  appropriate  choice,  if  the  remainder  of  the  rocket  can  reach  mission  objectives.  It  is 
possible  to  significantly  reduce  the  overall  length  of  the  aerospike  nozzle  while  avoiding 
significant  impact  to  the  overall  performance. 

Cooling  Jacket  Theory 

The  rocket  works  by  transferring  chemical  energy  into  useful  propulsive  energy  in 
the  combustion  chamber.  To  reach  this  condition,  the  pressure  and  temperature  of  the 
propellants  must  be  correct  for  efficient  combustion.  In  the  DEAN  concept,  a 
regenerative  cooling  expander  cycle  removes  energy  from  the  combustion  chamber  and 
nozzle  and  uses  that  energy  to  transport  and  pressurize  the  propellants  from  the  storage 
tanks  to  the  injector  and  chamber.  The  cooling  function  prevents  hardware  failure  due  to 
weakening  from  high  temperature  and  transfers  that  energy  to  the  working  fluids  via  the 
turbomachinery.  It  is  important  to  understand  how  the  energy  moves  through  the  system 
so  that  realistic  modeling  of  the  perfonnance  of  the  rocket  can  be  understood.  Energy  is 
conducted  through  the  walls  of  the  combustion  chamber  to  the  cooling  jacket  and  then  is 
imparted  to  the  heat  exchanger  working  fluids,  the  propellants,  via  convection. 

The  first  method  of  transfer  of  thermal  energy  is  convection  from  the  hot 
combustion  gases  to  the  walls  of  the  combustion  chamber  and  nozzle.  The  basic  equation 
to  model  this  heat  transfer  is  [7]: 

q  =  hg(Taw  -  Twg)  (20) 

where 

q  =  Heat  flux  (Btu/in2) 
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2 

hg  =  Gas-side  heat  transfer  coefficient  (Btu/(in"R)) 
Taw  =  Adiabatic  wall  temperature  of  the  gas  (R) 

Twg  =  Hot  gas  side  local  chamber  wall  temperature  (R) 


The  challenge  in  this  development  then  becomes  the  estimation  of  the  heat 
transfer  coefficient.  The  empirical  relation  developed  by  Bartz  is  used  to  estimate  hg. 
This  estimate  is  one  of  the  largest  sources  of  error  for  the  rocket  model  development  and 
is  on  the  order  of  ten  to  twenty  percent  [7]. 
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where 


Dt  =  Throat  diameter  (in) 

Cp  =  Constant  pressure  specific  heat  (Btu/lb  °F) 

Pr  =  Prandtl  number  =  pCp/k 
p  =  viscosity  (lb/in  s) 

R  =  Nozzle  radius  of  curvature  at  throat  (in) 
a  =  Correction  factor  for  property  variations  over  boundary  layer 
A  =  Area  along  chamber  axis  (in") 

Similarly,  the  thermal  energy  is  then  imparted  to  the  propellants  through 
convection  through  the  cooling  jacket  walls.  Again,  the  challenge  with  the  problem  is 
estimating  the  heat  transfer  coefficient,  now  hc  for  the  cold  side.  Concern  must  be  given 
to  the  phase  of  the  coolant  liquids.  If  the  pressure  of  the  fluid  is  below  the  critical 
pressure,  limited  vaporization  can  occur  as  nucleate  boiling,  that  is,  a  vapor  bubble  will 
form  in  the  liquid  and  grow  until  it  is  constrained  by  the  rate  of  condensation  at  the  edge 
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of  the  bubble.  The  heat  flux  into  the  cooling  liquid  can  increase  tremendously  with  a 
relatively  small  increase  in  the  coolant  side  wall  temperature  over  the  nucleate  boiling 
regime,  but  there  is  an  upper  limit  of  heat  flux  at  which  too  many  bubbles  will  form  and 
cause  the  coolant  side  wall  temperature  to  dramatically  increase,  usually  causing 
structural  failure.  The  upper  boundary  of  nucleate  boiling  can  therefore  be  used  as  a 
practical  design  boundary. 

If  the  cooling  fluid  is  above  the  critical  pressure,  no  boiling  can  occur.  Although 
the  high  rate  of  heat  flux  associated  with  nucleate  boiling  is  attractive,  the  initial  design 
choice  of  an  expander  cycle  indicates  supercritical  pressures  are  needed  to  avoid 
condensation  in  the  turbomachinery.  With  the  operating  regime  selected  as  supercritical 
pressures,  the  empirical  relation,  the  Sieder-Tate  equation  for  turbulent  flow  can  be  used 
to  approximate  hc. 


0.029 Cpii0  2  (G°-8\  ,Tco\ 

pr2/3  yd0.2  )  [TwJ 


(22) 


where 

2 

G  =  Coolant  weight  flowrate  per  unit  area  (lb/in  s) 
d  =  Coolant  passage  hydraulic  diameter  (in) 

Tco  =  Coolant  bulk  temp  (R) 

Twc  =  Coolant  side  wall  temperature  (R) 


The  total  cooling  capacity  of  the  regenerative  cooling  system  can  be  estimated  as: 

Qc  =  wcCp(Tcc  -  Tci)  (23) 

where 
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Qc  =  Coolant  capacity  (Btu/s) 


wc  =  Coolant  mass  flow  rate  (lb/s) 

Tcc  =  Coolant  critical  temperature  (R) 

Tcj  =  Coolant  inlet  temperature  (R) 
Oxygen  and  methane  critical  point  data  is  presented  in  Table  1. 


Table  1:  Propellant  Critical  Point  Data. 


Propellant 

Critical  Temperature  (R) 

Critical  Pressure  (psia) 

Oxygen 

278.26 

731.425 

Methane 

343.62 

672.927 

Hydrogen 

59.76 

188.108 

The  total  amount  of  energy  removed  from  the  chamber  and  imparted  into  the 
propellants  is  an  important  performance  parameter  where  hydrogen  has  a  major 
advantage  over  methane  due  to  heat  transfer  characteristics.  In  the  updated  methane 
DEAN  model,  it  will  be  necessary  to  verify  that  the  modeled  heat  transfer  rate,  Q  is 
safely  less  than  Qc  as  calculated  above. 

The  transfer  of  thermal  energy  from  the  inner  wall  of  the  combustion  chamber  to 
the  inner  wall  of  the  cooling  jacket  involves  conduction.  Conductive  heat  transfer  is 
defined  by  Fourier’s  Law  [4]: 


(24) 


where 


Q  =  Rate  of  heat  transfer  (Btu/s) 

K  =  Thennal  conductivity  (Btu/(ft*s*R)) 

A  =  Cross  sectional  area  perpendicular  to  flow  (in  ) 
T  =  Temperature  (R) 
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x  =  Distance  in  direction  of  heat  flow  (in) 


At  every  step  (chamber  to  wall,  through  the  walls,  wall  to  coolant),  the  rate  of 
heat  transfer  is  the  same,  so  it  is  possible  to  calculate  the  wall  temperatures  based  on 
material  properties.  The  desirable  characteristics  of  the  material  for  the  wall  are  high 
thermal  conductivity  and  a  high  maximum  possible  operating  temperature.  Additionally, 
high  pressures  found  in  the  cooling  channels  and  the  chamber  drive  structural 
requirements  of  the  hardware,  dictating  the  thickness  of  walls. 

Small  passages  with  high  flow  rates  and  high  pressures  lead  to  large  pressure 
drops  which  are  a  drain  on  the  efficiency  of  the  design  and  should  therefore  be 
minimized.  Therefore,  the  pressure  drop  across  the  channels  is  dependent  on  plumbing 
layout.  Changes  in  direction  and  contractions  and  expansions  must  be  minimized.  The 
pressure  drop  can  be  calculated  by: 
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where 

Ap  =  Pressure  drop  (lb/in") 

L  =  Length  of  the  portion  being  measured  (in) 
d  =  Equivalent  average  diameter  (in) 
p  =  Average  coolant  density  (lb/in  ) 

Vco  =  Coolant  flow  velocity  (in/s) 
g  =  Mass  conversion  factor  (in/s~) 
f  =  Friction  loss  coefficient 


26 


Materials 


As  noted  earlier,  the  combustion  process  creates  extreme  environments  which  can 
cause  failures  through  multiple  mechanisms.  The  different  elements  of  the  engine  must 
survive  the  associated  stresses  with  an  appropriate  safety  margin  (1.5  for  the  DEAN  [8]). 
High  temperatures  can  cause  structural  weakness  and  the  combustion  chamber  and  nozzle 
must  be  able  to  withstand  high  pressures  while  held  at  high  temperatures.  The  cooling 
channels  must  survive  highly  pressurized,  reactive  fluid  streams.  The  turbomachinery  is 
exposed  to  the  corrosive  fluids  and  thennal  gradients.  To  attain  the  reusability  and 
reliability  required  by  the  NGE  proposal,  all  structures  must  remain  below  material  limits 
with  margins  of  safety.  Martin’s  material  choices  were  chosen  based  on  processes  only 
reaching  50%  of  melting  points  [13].  Additionally,  Martin  accounted  for  material 
compatibility  with  the  propellants  of  liquid  hydrogen  and  liquid  oxygen.  The  cooled 
section  of  the  aerospike  was  copper  and  the  tip  was  made  of  niobium.  Silicon  carbide 
comprised  the  cooled  section  of  the  chamber  while  the  jacket  was  made  of  aluminum. 

Upon  more  detailed  investigation  into  the  optimization  of  thrust  to  weight,  Hall 
selected  materials  to  meet  the  NGE  goals  and  ease  manufacturability  [8],  Oxygen  free 
copper  was  selected  for  all  of  the  aerospike  components.  The  chamber  again  used  silicon 
carbide  because  of  the  favorable  thennal  properties.  To  avoid  reaction  with  the 
propellants,  the  hydrogen  plumbing  was  made  from  INCOLOY  909  and  the  oxygen 
plumbing  was  made  from  INCONEL  718. 
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Previous  Work 


Martin 

David  Martin  used  NASA's  NPSS  to  design  the  first  DEAN  model,  which  is  built 
up  through  different  elements,  each  representing  the  different  components  of  the  rocket 
design  [13].  Martin's  methodology  assumed  rocket  engine  parameters  based  on 
historically  successful  designs  and  used  NPSS  to  perfonn  a  power  balance  and  a  pressure 
budget.  NPSS  validated  the  reasonableness  of  the  initial  estimates  and  calculated  more 
accurate  parameters.  Once  the  model  solution  closed,  Martin  explored  design  changes  to 
increase  performance. 

Martin's  objectives  were  to  explore  the  feasibility  of  meeting  the  NGE  research 
goals  with  the  hydrogen  DEAN  concept,  establish  a  design  process  focused  on  the  energy 
conversion  sections  of  the  rocket,  and  to  perform  detailed  design  analysis  of  energy 
transfer  components.  The  results  in  Table  2  demonstrate  that  the  DEAN  design  was  an 
excellent  candidate  for  further  study  in  pursuit  of  the  NGE  goals. 


Table  2:  Martin’s  DEAN  Design 


DEAN 

DEAN  Goals 

Comparison 

Vacuum  Thrust  (lbf) 

57,231 

50,000 

+  14.5% 

Vacuum  Isp  (s) 

472 

464 

+  1.7% 

Additionally,  Sierra  Engineering's  software  package  Two  Dimensional  Kinetics 
(TDK’04™)  was  used  to  develop  the  aerospike  geometry.  TDK’04  used  combustion 
chamber  data  (pressure  and  temperature)  to  develop  the  nozzle  contour.  Further  work 
explored  performance  loss  from  aerospike  truncation,  with  the  goal  of  volume  and  weight 
savings.  The  result  of  this  analysis  was  that  significant  truncation  of  the  nozzle  (up  to 
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50%)  past  the  throat  resulted  in  minor  decreases  (on  the  order  of  2%)  in  system 
performance. 

Simmons 

Joseph  Simmons  upgraded  the  DEAN  model  with  a  focus  on  automation  of 
parametric  studies  [14].  The  goal  of  the  work  was  to  develop  the  tools  necessary  to 
optimize  the  thrust  to  weight  ratio  of  the  DEAN  model,  and  this  was  accomplished 
through  adding  geometry  parameters  into  the  model  and  exploring  the  O/F  ratio.  The 
result  was  a  more  robust  DEAN  modeling  tool  that  produced  a  design  with  similar 
performance,  but  a  significantly  shorter,  and  therefore  lighter,  engine.  To  accomplish 
this,  the  NPSS  independent  variables  and  dependent  conditions  were  adjusted.  The  NPSS 
model  was  wrapped  in  ModelCenter  for  greater  control  and  automation.  The  input  and 
output  variables  of  NPSS  were  therefore  readily  available  for  study  and  additional 
analysis. 

Simmons'  work  extracted  top  level  parameters  from  the  inner  workings  of  the 
DEAN  model  to  make  them  available  for  parametric  evaluation.  For  example,  before 
Simmons'  work,  the  O/F  ratio  was  controlled  indirectly  through  manipulating  pump 
pressure  ratios.  By  reallocating  the  O/F  ratio  as  a  design  variable,  Simmons  increased 
performance  by  decreasing  the  O/F  ratio  from  7  to  6,  thereby  finding  an  optimum  value 
for  this  engine  and  cycle.  An  optimum  O/F  ratio  increases  the  amount  of  energy  that  can 
be  harnessed  from  the  combustion  reaction  and  therefore  may  yield  increased 
perfonnance,  provided  chamber  and  throat  temperatures  stay  within  material  temperature 
limits.  The  expander  cycle  active  cooling  of  the  surfaces  in  question  should  allow  for  this 
increase  in  perfonnance,  and  the  DEAN  model  checks  for  these  limits  through  energy 
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balances  and  material  property  checks.  The  ModelCenter  controlled  NPSS  model 
produces  a  closed  design  for  each  set  of  input  parameters  so  geometries,  chamber 
pressures,  etc.  may  all  change  between  designs;  for  this  reason  it  was  important  to 
reorganize  the  input  and  output  variables. 

The  end  result  of  the  effort  was  the  development  of  a  more  accurate  model  with 
an  increase  in  the  number  of  customizable  user  inputs.  The  automation  of  the  parametric 
analyses  provides  the  user  with  the  ability  to  quickly  create  a  family  of  rocket  engines  at 
various  discreet,  closed  designs.  The  generated  designs  can  then  be  analyzed  for  desired 
characteristics.  Using  Martin's  design  input  characteristics,  the  new  design  process 
developed  an  engine  with  the  following  performance. 


Table  3:  Simmons’  Parametric  Results. 


Initial  DEAN 

DEAN  -  v2 

Requirements 

Vacuum  Thrust  (lbf) 

57,231 

50,000 

50,000 

Vacuum  Isp  (s) 

472  s 

464 

464 

O/F 

7 

6 

N/A 

Mass  Flow  (lbm/s) 

121 

104 

N/A 

Length  (in) 

37.2 

27.9 

N/A 

Hall 

Hall  used  the  updated  ModelCenter/NPSS  DEAN  design  tool  to  develop  a 
substantially  more  detailed  system  level  design  point  through  the  use  of  automated  trades 
and  parametric  studies  [8],  Additionally,  the  reusability  of  the  design  was  examined 
through  a  materials  survey  which  provided  material  property  data  to  the  model.  The 
temperature  dependent  material  property  data  was  incorporated  into  automated  checks  for 
structural  integrity.  Different  sections  of  the  engine  experience  different  stresses  and 
temperatures.  The  performance  calculations  Hall  developed  estimate  the  volumes  of  the 
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rocket  components  (and  therefore  mass)  and  check  the  resulting  pieces  against  selectable 
failure  criteria.  The  resulting  design  perfonnance  is  listed  in  Table  4.  Hall’s  additions  to 
the  design  tool  included  much  more  accurate  modeling  of  the  aerospike  nozzle  behavior, 
including  losses  and  reasonable  expansion  ratios.  Additionally,  Hall  finalized  the 
capability  to  analyze  vacuum  thrust-to-weight  ratio,  giving  the  modeling  tools 
substantially  more  utility  in  the  pursuit  of  the  NGE  requirements. 


Table  4:  Hall's  DEAN  Design  Performance. 


DEAN  -  v3 

DEAN  Goals 

Comparison 

Vacuum  Thrust  (lbf) 

50,161 

50,000 

+0.3% 

Vacuum  T/W 

142.2 

106.5 

+33.5% 

Vacuum  Isp  (s) 

430.6 

464 

-7.3% 

Preliminary  Design  Decisions 

Dual  Expander  Cycle 

Traditional  single  expander  rocket  systems  use  fuel  flow,  usually  liquid  hydrogen, 
to  remove  heat  from  the  combustion  chamber  and  nozzle.  The  heated  fuel  then  flows 
through  a  turbine  to  provide  mechanical  work  to  the  fuel  and  oxidizer  pumps.  The  pumps 
pressurize  the  fuel  and  oxidizer  for  the  combustion  chamber  injector  system.  The  dual 
expander  system  differs  in  that  both  the  fuel  and  oxidizer  are  used  as  working  fluids; 
therefore  the  design  includes  a  fuel  expander  cycle  and  a  completely  separate  oxidizer 
expander  cycle.  By  separating  the  two  flow  systems,  the  likelihood  of  failure  due  to 
fuel/oxidizer  mixing  prior  to  the  combustion  chamber  is  reduced.  Also,  the  separation 
increases  options  for  material  selection  for  either  feed  system  based  on  the  material 
compatibility  of  the  fluid.  In  addition,  by  using  essentially  the  full  mass  flow  of  the 
engine  to  power  the  pumps,  the  pump  pressure  ratio  can  be  increased  and  the  turbine 
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pressure  ratio  decreased,  permitting  increased  chamber  pressure.  Increased  chamber 
pressure  for  a  given  propellant  flow  rate  increases  performance  by  lowering  the  area  of 
the  throat  and  all  subsequently  derived  engine  dimensions,  cutting  engine  mass  [4]. 

Aerospike  Nozzle 

The  first  advantage  of  the  aerospike  nozzle  is  geometry  that  is  conducive  to  a  dual 
expander  design.  Because  the  combustion  flow  occurs  in  an  annular  region,  there  are  two 
separated  surfaces  available  for  heat  transfer  from  the  chamber  into  the  pre-combustion 
working  fluids.  This  separation  should  ease  plumbing  design  and  therefore  lower  total 
plumbing/fuel  system  weight.  Savings  would  be  realized  in  a  specific  stage  design, 
which  is  beyond  the  scope  of  this  project.  Second,  the  aerospike  nozzle  can  operate 
optimally  at  all  altitudes  of  flight.  Although  this  benefit  seems  of  little  consequence  for 
an  upper  stage  engine,  developmental  test  and  evaluation  would  be  simplified.  Because 
of  the  high  cost  of  space  systems  acquisitions,  space  launch  customers  demand  a 
rigorously  tested  launch  solution.  This  design  choice  therefore  has  the  potential  to 
dramatically  reduce  the  cost  of  test  and  evaluation  of  a  launch  acquisition  by  reducing  the 
complexity  of  the  required  testing  infrastructure,  as  all  developmental  testing  could  be 
done  at  local  atmospheric  conditions,  requiring  a  low  pressure  chamber  only  for  final 
verification  testing. 

Methane  Fuel 

Some  properties  of  cryogenic  methane  are  appealing  when  compared  to  hydrogen. 
The  significantly  higher  boiling  point  temperature  and  the  larger  enthalpy  of  vaporization 
at  the  normal  boiling  point  (NBP)  imply  that  there  could  be  a  gain  of  efficiency  on  the 
pad  waiting  for  launch  (See  Table  5).  Hydrogen  boil-off  on  the  pad  is  estimated  at  1.2% 
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per  hour,  requiring  crews  to  “top  off’  launch  vehicles  near  launch  [15].  Less  energy  will 
be  needed  to  keep  the  cryogenic  fuel  at  an  appropriate  temperature,  the  rate  of  energy 
conduction  into  insulated  tanks  will  be  reduced  (because  of  the  smaller  change  in 
temperature),  and  the  greater  enthalpy  of  vaporization  means  there  should  be  less  boil  off 
on  the  pad,  reducing  the  extra  fuel  to  be  added  in  a  design  through  engineering  margin. 


Table  5:  Properties  of  Cryogenic  Fuels  [16]. 


Property 

Unit 

Hydrogen 

Methane 

Molecular  Weight 

g/mole 

2.0159 

16.043 

Tb  @  1  atm  or  NBP 

K 

20.28 

111.668 

Enthalpy  of  Vap  @  NBP 

J/R 

445 

510.83 

Liquid  Density  @  NBP 

g/mL 

0.0708 

0.4224 

Heat  Capacity  (1)  @  NBP 

J/gK 

9.668 

3.481 

The  combustion  products  of  methane  have  a  relatively  low  molecular  weight.  The 
chemical  species  of  greatest  abundance  in  methane/oxygen  combustion  products  are  H2O, 
CO,  H2  ,  and  CO2  .  The  relative  amount  of  product  species  depends  on  several  factors  to 
include  combustion  chamber  geometry,  temperature,  and  pressure  [12].  In  equations  26 
through  28,  it  is  shown  that  Isp  is  inversely  related  to  products  molecular  weight,  MW  [4]. 
While  methane  compares  unfavorably  with  hydrogen,  it  compares  favorably  against 
many  other  fuels. 
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Toxicity  of  combustion  products  is  not  terribly  important  for  upper  stage  engines, 
but  low  toxicity  lends  itself  to  the  development  of  a  reusable,  safe,  and  very  testable 
device.  When  combined  with  the  atmospheric  compensation  of  the  aerospike  nozzle,  the 
DEAN  concept  is  relatively  easy  to  test  on  the  ground. 


Basis  of  Comparison 

Considerable  work  has  been  applied  to  the  development  of  aerospike  engines  in 
the  past,  so  it  is  important  to  note  the  perfonnance  of  previous  work  to  detennine  if 
current  efforts  provide  a  comparative  advantage.  The  first  opportunity  for  comparison  is 
the  previous  work  of  Hall,  a  similarly  derived  DEAN  with  hydrogen  fuel.  Unfortunately, 
this  will  not  quite  provide  a  perfect  comparison  as  Hall’s  requirements  were  somewhat 
different  than  what  is  currently  being  investigated.  The  different  fuel  will  give  the 
hydrogen  based  design  a  higher  Isp  and  total  performance  advantage,  but  the  target  of 
50,000  lbf  [8]  of  thrust  prevents  a  direct  comparison. 

Previous  methane  based  engine  concepts  exist  in  the  literature.  Klepikov  et  al. 
present  a  variety  of  engines  employing  methane  as  fuel  [6].  The  first,  RD-185  is  an  upper 
stage  engine  designed  for  use  with  methane.  The  second  two,  RD-167  and  RD-160,  are 
modifications  of  the  kerosene  engines,  RD-134  and  RD-161,  respectfully.  All  three  of 
these  engines  include  bell  nozzles  with  large  area  ratios.  The  perfonnance  parameters  of 
several  historical  engines  of  interest  are  in  Table  6. 
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Table  6:  Comparable  Engine  Performance. 


Engine 

RD-185[6] 

RD-167161 

RD-160[6] 

NGE  -  Derived 

nge[5] 

H-DEAN181 

Fuel 

ch4 

ch4 

ch4 

ch4 

h2 

h2 

Thrust  (lbf) 

40,344 

79,366 

4,409 

25,000 

25,000 

44,694 

Up  (s) 

378 

379 

380.6 

383 

465 

429.8 

T/W 

44 

63 

15.5 

108 

N/A 

142.2 
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III.  Methodology 


With  research  goals  defined,  this  chapter  describes  the  method  by  which  the 
MDEAN  model  was  developed.  The  chapter  begins  with  a  discussion  of  the  different 
software  tools  used  in  the  development  and  includes  descriptions  of  some  of  the  critical 
software  functions  necessary  for  model  calculations.  Next  is  a  more  focused  description 
of  how  NPSS  uses  thermochemistry  and  an  explanation  on  the  development  of  the 
different  fluid  property  tables.  Finally,  there  is  a  discussion  on  the  modifications  to 
NPSS  and  ModelCenter  that  were  necessary  for  the  hydrogen  to  methane  fuel  conversion. 

Tools 

ModelCenter 

ModelCenter  by  Phoenix  Integration  is  a  program  designed  for  the  development 
and  study  of  software  models.  Many  other  programs  can  be  “wrapped”  by  ModelCenter 
to  control  input  and  to  accept  and  analyze  output.  A  key  feature  is  the  automation  of 
parametric  studies,  allowing  for  a  large  amount  of  untended  software  execution  that  is 
necessary  to  explore  a  design  space.  Additionally,  ModelCenter  has  a  design-of- 
experiments  function  that  allows  the  user  to  statistically  analyze  the  response  of  a  system 
for  efficient  optimization  of  input  design  variables.  The  visualization  suite  provides  a 
convenient  method  of  displaying  and  analyzing  data. 

ModelCenter  is  the  top  level  software  tool  used  to  control  the  command  line 
driven  modeling  programs,  but  it  also  has  several  features  that  have  aided  in  the 
development  and  execution  of  the  DEAN.  A  deeper  level  of  automation  and  control  of 
the  NPSS  model  can  be  realized  through  the  use  of  ModelCenter’ s  included  scripting 
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tools.  Both  VBScript  and  JavaScript  were  used  in  the  DEAN  models  to  perform 
intennediate  calculations  and  provide  functionality  checks  on  the  model.  Additionally, 
the  data  import  function  allows  a  user  to  import  data  sets  for  analysis.  The  RMS  toolkit 
can  then  be  used  to  provide  a  polynomial  or  Kriging  regression  of  data  sets.  The 
MDEAN  primarily  uses  Kriging  regression  functions  for  thennodynamic  tables  based  on 
enthalpy,  which  will  be  discussed  in  detail  later. 

Kriging  Estimator 

The  Kriging  function  is  a  method  of  interpolation  which  relies  on  the  observation 
of  surrounding  data  points  weighted  for  spatial  covariance  [17].  This  method  of  the 
linear  regression  estimator,  Z*(u)  is  defined: 

n(u ) 

Z*{u)  -  m(u )  =  ^  Aa  [Z (it)  -  m(ua )]  (29) 

OC=l 

where 

u,  ua  =  location  vector  for  estimation  point  and  neighboring  points  with  index  a 
n(u)  =  number  of  neighborhood  points  used  for  Z*(u)  estimation 
m(u),  m(ua)  =  expected  mean  values  of  Z(u)  and  Z(ua) 
ka(u)  =  Kriging  weight;  each  estimated  Z(ua)  will  have  a  different  weight 
The  Kriging  weight,  A„,  is  the  key  to  this  method  of  estimation  and  is  assigned  to 
neighboring  data  points  to  find  an  estimate  at  the  queried  point  u.  The  derivation  to 
calculate  the  weights  and  an  example  can  be  found  in  reference  [17].  The  weight  of  a 
point  often  decreases  as  the  distance  from  the  location  (temperature  and  pressure  for  the 
purposes  of  the  MDEAN)  to  be  estimated  and  that  point  increase.  The  ModelCenter 
Kriging  function  handles  all  calculations  for  the  user,  but  care  must  be  taken  in  the  data 
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selection  used  in  the  Kriging  generation.  As  with  all  interpolation  algorithms,  evenly 
spaced  data  covering  all  of  the  dataset  to  be  analyzed  will  provide  the  best  results.  Step 
sizes  must  be  small  enough  to  provide  ample  coverage  of  non-linear  behavior  (or  simply 
fast  changing  behavior),  but  the  ModelCenter  algorithm  indicates  best  results  occur  with 
fewer  than  one  thousand  data  points,  so  there  is  a  tradeoff  between  step  size  and 
accuracy.  Because  the  Kriging  function  is  a  weighted  average,  it  will  often  produce 
estimates  that  are  very  smooth,  sometimes  more  so  than  the  actual  data  [17].  Another 
characteristic  of  this  algorithm  is  that  it  should  provide  exact  results  if  a  data  point  used 
in  Kriging  generation  is  then  queried  with  the  function.  This  makes  sense  as  the  data 
point  would  be  given  a  weighting  to  the  exclusion  of  other  points.  Although  it  is  positive 
that  the  Kriging  can  perfectly  estimate  some  discreet  data  points,  this  phenomenon 
increases  the  workload  with  regards  to  error  analysis  as  more  source  data  is  required  to 
check  the  error  of  the  Kriging  generated  estimates. 

Polynomial  Regression 

In  addition  to  Kriging  interpolation,  ModelCenter  provides  a  polynomial 
regression  tool  to  develop  models  for  data  sets.  The  polynomial  regression  is  another 
method  of  linear  regression,  but  the  result  is  a  polynomial  function  that  approximates  the 
shape  of  the  input  data.  While  the  fit  of  this  regression  was  usually  lower  than  the 
Kriging,  as  evidenced  by  R  values,  the  method  proved  useful  for  developing  one  three 
dimensional  fluid  property  table.  The  polynomial  function  is  advantageous  because  it 
can  be  used  directly  in  the  place  of  the  interpolation  table.  This  is  convenient  as  it 
eliminates  the  error  from  the  linear  interpolation.  Unlike  the  Kriging,  the  polynomial 
function  doesn’t  exactly  predict  the  input  data,  so  error  analysis  can  largely  be  done  with 
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the  input  data.  Also,  this  linear  regression  can  accommodate  significantly  more  data,  so 
error  can  be  reduced  in  areas  of  interest  by  adding  input  data. 

Numerical  Propulsion  Simulation  System  (NPSS) 

NPSS  is  a  NASA  developed  simulation  program  used  to  model  engines  and 
provide  estimates  of  performance  and  is  at  the  heart  of  the  DEAN  model.  While  NPSS 
was  originally  developed  for  air  breathing  applications,  rocket  based  elements  have  been 
developed  as  of  release  1.65  and  the  iterative  solution  methodology  of  the  program  lends 
itself  to  solving  rocket  systems. 

Elements 

NPSS  works  by  combining  elements  that  are  mathematical  representations  of  unit 
operations  or  physical  phenomena.  The  NPSS  Rockets  Supplement  [18]  lists  the  relevant 
input  and  output  variable  names  and  units.  Each  element  is  an  independent  file  of 
software  which  can  be  found  in  the  InterpComponents  directory  under  the  NPSS 
installation  directory  (for  version  1.65).  Examining  this  code  is  helpful  for  understanding 
the  underlying  equations  that  are  calculated  by  the  solver.  Some  elements,  such  as  the 
combustion  chamber,  have  independent  variables  and  dependent  conditions  associated 
with  them  which  will  be  utilized  for  the  overall  system  solver  if  the  solver  auto  setup 
function  is  employed.  Independent  variables  and  dependent  conditions  are  discussed 
below. 

Ports 

Ports  are  the  software  method  by  which  elements  are  linked  and  include  fluid 
ports,  heat  flow  ports,  or  mechanical  energy  ports.  Each  element  will  have  ports  that 
must  be  connected  to  the  ports  of  the  other  elements  to  provide  the  flow  of  mass  or 
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energy.  The  linkPorts  function  builds  the  network  of  connections  required  by  the 
iterative  solver.  A  port  is  a  software  linkage  and  does  not  represent  a  physical 
connection.  A  fluid  port  does  not  include  an  associated  pressure  drop  due  to  frictional 
losses.  It  is  therefore  necessary  to  include  an  element  that  has  a  frictional  pressure  loss  to 
represent  the  physical  fluid  connection.  That  element  (Valve04  in  DEAN)  will  have  at 
least  two  ports  which  must  be  linked  to  other  elements. 

Solver 

The  NPSS  solver  is  an  iterative  algorithm  which  attempts  to  satisfy  dependent 
conditions  through  the  adjustment  of  independent  variables.  The  goal  of  the  solver  is  to 
drive  the  system  model  to  a  consistent,  converged  state.  Several  elements  have  built  in 
independents/dependents  and  are  usually  associated  with  the  fundamental  continuity 
equations  that  are  to  be  expected  of  a  flow  system.  For  example,  mass  flow  entering  an 
element  is  equal  to  mass  flow  exiting  the  same  element.  It  is  inevitable  a  model  will 
require  a  controlling  variable  or  condition  that  is  not  inherent  in  the  comprising  elements, 
so  a  user  can  add  independent  variables  and  dependent  conditions.  Independent  variables 
and  dependent  conditions  are  associated  with  the  system’s  solution  method  and  therefore 
can  only  be  added  in  pairs,  unless  the  existing  variables/conditions  are  modified,  i.e.  there 
need  to  be  an  equal  number  of  independents  and  dependents  so  the  solver  can  work  with 
a  square  matrix. 

The  dependent  conditions  are  organized  as  an  equation,  left  hand  side  =  right  hand 
side  (in  the  example  above,  mass  flow  in  =  mass  flow  out).  Each  side  of  this  equation  is 
dependent  on  one  or  more  of  the  independent  variables.  For  the  sake  of  this 
development,  consider  each  side  being  dependent  on  the  independent  variable  x.  The 
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inequality  defines  an  error  term  which  is  equal  to  the  difference  of  the  two  sides  of  the 
equation,  but  this  simple  error  term  is  undesirable  as  the  tolerance  can  then  be  variable 
amongst  the  different  equations/dependents.  For  example,  if  each  side  of  the  dependent 
equation  drives  to  a  value  of  1,  a  tolerance  of  10'5  might  be  appropriate,  but  if  each  side 
of  the  dependent  equation  is  driving  to  10,000,  then  the  same  tolerance  might  be  too 
restrictive  [19].  Therefore,  the  error  tenn  is  equal  to  the  difference  of  the  two  sides  of  the 
equation  divided  by  a  reference  term,  which  should  be  the  same  order  of  magnitude  as  the 
calculated  solution  to  the  equation.  A  quasi-Newton  method  is  then  used  to  drive  the 
error  term  to  zero.  This  means  that  the  error  term  is  calculated  at  the  guessed  value  of  the 
independent  variable,  x,  and  then  the  derivative  of  the  error  tenn  curve  is  calculated  at  the 
guessed  value  of  interest  of  x.  This  derivative  at  x  forms  a  line  that  is  propagated  to 
where  the  enor  term  is  zero,  and  the  value  of  x  is  found  at  this  point.  The  process  begins 
again  along  the  error  curve  for  the  new  value  of  x.  The  quasi  nature  of  this  method 
means  the  true  slope  is  not  calculated  for  every  iteration  in  an  effort  to  save  computation 
time.  This  process  is  conducted  amongst  all  independent/dependent  interactions  via  a 
matrix  known  as  the  Jacobian.  The  Jacobian  matrix  is  updated  after  a  number  of 
iterations  via  Broyden’s  method,  which  is  explained  in  the  NPSS  User’s  Guide  [19]. 

Output 

As  stated  before,  NPSS  uses  C++  syntax  and  therefore  uses  the  cout  command  to 
print  data.  The  cout  function  is  used  extensively  in  the  MDEAN  model  to  output  the 
required  NPSS  output  variables  to  a  text  file  for  review  by  a  user  or  ModelCenter.  A 
limitation  of  the  cout  command  with  regards  to  NPSS  is  it  only  prints  converged, 
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successful  runs  or  the  final  state  before  run  failure.  A  different  data  dump  function  would 
be  useful  to  see  the  evolution  of  an  iterating  run. 

Chemical  Equilibrium  with  Applications  (CEA) 

CEA  is  a  NASA  developed  program  used  to  calculate  thermodynamic  properties 
at  different  states.  The  primary  use  of  CEA  for  the  current  work  is  for  the  calculation  of 
thermochemical  properties  of  the  combustion  products  [20],  The  CEA  rocket  problem 
function  has  several  options  which  are  set  to  specify  the  rocket  parameters  most  desirable 
for  the  modeling  application.  For  example,  the  DEAN  combustion  products  are  modeled 
by  specifying  the  chamber  state  (pressure  and  temperature),  the  type  of  equilibrium 
(frozen  flow),  the  O/F  ratio,  and  the  reactants  (methane  and  oxygen).  CEA  assumes  an 
ideal  gas  for  the  equation  of  state  of  the  mixture  of  chemical  species  found  in  the  reaction 
products.  Chemical  equilibrium  is  then  defined  by  the  minimization  of  free  energy. 

The  following  mathematical  derivation  is  from  the  CEA  guide,  NASA  Reference 
Publication  1311  [21].  Gibb’s  energy  per  kilogram  mixture  for  NS  species  is  given  by 
equation  30  where  n,  is  the  number  of  kilogram-moles  of  species  j  per  kilogram  of 
mixture: 


9 


NS 

= 


(30) 


j= i 


And  iij,  the  chemical  potential  per  kilogram-mole  of  species  j  is: 


9j 


'dg\ 

^dUj  )  T,P,ni+J 


(31) 


The  equilibrium  condition  is  described  by  the  constrained  minimization  of  the 
Gibb’s  energy,  the  constraint  being  a  mass  balance  for  the  element  i,  or 
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(32) 


NS 

dijiij  —  bf  —  0  (i  =  1, ... ,  Z) 

7  =  1 

where  the  first  tenn  is  the  number  of  kilogram-atoms  of  element  i  per  kilogram  of 
mixture  (the  number  of  element  i  atoms  in  molecule  j)  and  b  is  the  assigned  number  of 
kilogram-atoms  of  element  i  per  kilogram  of  total  reactants  (total  number  of  atoms  of 
element  i ).  The  index  /  is  the  number  of  chemical  elements.  This  constrained 
minimization  problem  is  solved  via  Lagrange  multipliers  with  the  following  definition: 

l  ( NS  \ 

G  =  g  +  ^  h  (  ^  dijUj  -  bf  J  (33) 

i= i  \J= 1  J 

are  the  Lagrange  multipliers.  The  method  of  Lagrange  multipliers  is  an 
optimization  technique  to  find  a  local  maxima  or  minima  subject  to  constraints,  the 
constraint  being  the  mass  balance.  The  equilibrium  conditions  are  then  calculated  by: 

NS  /  l  \  l 

5G  =  ^  Lj  +  ^  Aifly  j  Snj  +  -  bf)5Ai  =  0  (34) 

7=1  V  i=l  /  i=l 

And  as  the  variations  are  independents: 

i 

[ij  +  ^  Aidij  —  0  0  =  1-  - .  NS)  (35) 

i=l 

This  set  of  equations  can  then  be  solved  for  equilibrium  concentrations  with 
knowledge  of  the  chemical  potentials,  p,  [21]. 
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NPSS  Thermochemistry 


NPSS  supports  a  number  of  chemical  species  to  include  hydrogen  and  oxygen, 
but  methane  and  its  combustion  products  are  not  directly  supported.  There  are  several 
methods  of  incorporating  thennodynamic  data  into  the  DEAN  model,  and  each  different 
method  is  specified  in  the  model  with  the  setThermoPackage  command.  The  first  method 
explored  was  direct  integration  of  CEA  calls  into  the  NPSS  code.  When  a  calculation 
requires  thermodynamic  data,  NPSS  can  launch  CEA  and  acquire  the  necessary 
information.  Unfortunately,  the  employed  rocket  elements  were  developed  with 
unreacted  fluid  stations  and  therefore  cannot  make  use  of  the  CEA  thennodynamics 
package.  After  consulting  with  Wolverine  Ventures,  the  decision  was  made  to  develop 
custom  fluid  property  tables  for  the  chemical  species  of  interest.  A  fluid  property  table 
includes  several  thennodynamic  properties,  each  as  a  function  of  the  state  of  the  fluid  and 
can  be  seen  in  Table  7  [22]. 


Table  7:  Fluid  Property  Table  Supported  Parameters. 


Property 

Unit 

Property 

Unit 

Density 

lbm/ft3 

Viscosity 

lbm/(ft*sec) 

Enthalpy 

Btu/lbm 

Thennal  Conductivity 

Btu/(sec*ft*R) 

Internal  Energy 

Btu/lbm 

Ratio  of  Specific  Heats 

N/A 

Entropy 

Btu/(lbm*R) 

Prandtl  Number 

N/A 

Const  Vol  Heat  Cap 

Btu/(lbm*R) 

Molecular  Weight 

lbm/lb-mole 

Const  Press  Heat  Cap 

Btu/(lbm*R) 

Total  Gas  Constant 

Btu/(lbm*R) 

For  the  liquid  propellants,  the  National  Institute  of  Standards  and  Technology 
(NIST)  thermodynamic  database  was  used  as  the  source  of  data  [23].  Data  was  taken 
over  the  ranges  of  interest  for  methane  and  oxygen  with  ranges  and  step  sizes  detennined 
through  logical  analysis  of  the  problem  at  hand.  For  both  propellants,  the  phase  of  the 
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chemical  species  dictated  the  area  of  interest  for  thermodynamic  data.  Methane  and 
Oxygen  exist  in  the  rocket  as  pure  species  from  the  propellant  tanks  to  the  injector  and 
should  only  exist  as  liquids  or  supercritical  fluids  because  vapor  would  cause  flow 
instabilities  and  possible  cavitation  in  turbo-machinery.  The  phase  data  of  the  pure 
propellants  therefore  partially  dictate  the  ranges  of  temperatures  and  pressures  at  which 
data  should  be  tabulated.  Below,  Figure  5  shows  the  boiling  temperatures  of  methane 
over  a  range  of  pressures.  This  curve  represents  the  lower  boundary  of  the  query  for 
NIST  data.  Below  the  indicated  pressure,  the  pure  chemical  is  a  gas  and  is  therefore  not 
applicable  to  the  fluid  flow  problem. 


Figure  5:  Methane  Boiling  Point  Temp  vs.  Pressure  (Antoine  Equation). 


NIST  includes  very  high  pressure  data  for  both  oxygen  and  methane  with  upper 
limits  being  1 1,893  psia  and  145,037  psia  respectively.  For  the  purposes  of  the  DEAN, 
the  upper  boundaries  for  the  fluid  property  tables  are  significantly  reduced  based  on 
historical  data.  First,  the  upper  end  of  chamber  pressures  for  Hall’s  designs  is 
approximately  3000  psia  [8].  Accounting  for  an  approximate  30%  [4]  pressure  drop  for 
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an  injector  and  a  20%  drop  for  a  turbine  sets  a  reasonable  upper  pressure  bound  of  4600 
psia  for  all  propellants.  For  margin,  the  tables  were  extended  to  6000  psia  in  the 
supercritical  region.  If  a  fluid  station  has  a  state  beyond  the  limits  of  the  table,  NPSS  will 
extrapolate  a  value  for  a  queried  property.  This  represents  a  possible  source  of  error  and 
should  be  avoided  if  at  all  possible. 

The  temperature  limits  used  for  the  fluid  property  tables  are  outlined  in  Table  8. 
The  low  temperature  was  set  at  a  convenient  integer  slightly  above  the  species  melting 
point  at  one  atmosphere.  This  value  was  chosen  because  frozen  propellants  are 
undesirable  and  will  not  occur  in  the  DEAN  model  by  design.  The  NIST  database  upper 
temperature  limits  are  1 125  Rankin  for  methane  and  1800  R  for  oxygen  [23].  The  tables 
include  high  temperature  data  up  to  the  NIST  limit  in  an  attempt  to  eliminate  error  from 
extrapolation.  One  mode  by  which  the  NPSS  model  failed  to  converge  was  due  to  the 
oxygen  temperature  exceeding  the  upper  temperature  limit  of  the  fluid  property  table. 

The  extrapolation  method  was  set  to  Lagrange2  at  the  time  and  NPSS  returned  a  Cp  of 
approximately  -0.5.  This  caused  the  heat  transfer  coefficient  on  the  cooling  channel  side 
to  be  a  complex  number  (via  the  Bartz  estimate,  see  equation  21)  and  NPSS  was  unable 
to  continue  calculation.  The  temperature  range  of  the  oxygen  FPT  was  then  expanded  to 
the  extent  NIST  provides  and  the  extrapolation  method  was  changed  to  linear  for  all 
tables. 

Table  8:  Fluid  Property  Input  Data  Independent  Variable  Ranges. 


Methane 

Oxygen 

T,ow  (R) 

165 

100 

Thigh  (R) 

1115 

1800 

PioW  (psia) 

15 

15 

Phigh  (psia) 

6000 

5000 
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The  NPSS  User’s  Guide  [19]  dictates  the  fluid  property  table  syntax.  For  a  given 
temperature  (or  enthalpy  or  entropy),  a  pressure  sweep  and  the  corresponding  tabulated 
thermodynamic  properties  are  listed.  This  organization  dictates  attempting  to  record 
significant  pressure  sweeps  at  discreet  temperatures.  Fortunately,  the  NIST  database 
allows  the  setting  of  range  and  step  size  for  data  queries.  The  interpolation  method  for 
the  fluid  property  tables  is  linear.  This  limitation/feature  implies  a  requirement  for  the 
step  size  of  the  collected  data.  Over  the  ranges  of  the  fluid  properties  that  are  linear,  a 
large  step  size  can  be  used.  Over  ranges  where  nonlinear  behavior  is  observed  (near  the 
critical  point),  a  smaller  step  size  is  needed  to  reduce  interpolation  error. 

Methane  FPT  Generation 

Based  on  the  FPT  example  built  with  the  NPSS  Rocets  thermodynamic  package, 
the  methane  fluid  property  table  requires  the  following  individual  property  tables: 


Table  9:  Properties  Required  for  Methane  FPT. 


Property,  NPSS  Symbol 

Unit 

Independent  Variables 

Indeps  Call 

Density,  rho 

Ibn/fV 

h,P 

N/A 

Internal  Energy,  u 

Btu/lbm 

h,  P 

N/A 

Entropy,  s 

Btu/lbm/R 

h,  P 

N/A 

Viscosity,  mu 

lbm/ft/s 

h,P 

N/A 

Const  V  Specific  Heat,  Cv 

Btu/lbm/R 

h,P 

N/A 

Const  P  Specific  Heat,  Cp 

Btu/lbm/R 

h,P 

N/A 

Thermal  Conductivity,  k 

Btu/s/ft/R 

h,P 

N/A 

Ratio  of  Specific  Heats,  gam 

N/A 

h,P 

N/A 

Temperature,  T_h 

R 

h,P 

Thlndeps 

Enthalpy,  h_T 

Btu/lbm 

T,  P 

hTindeps 

Temperature,  T_s 

R 

s,P 

Tslndeps 

Enthalpy,  h_s 

Btu/lbm 

s,P 

hslndeps 

Total  Gas  Constant,  R 

Btu/lbm/R 

h,P 

N/A 
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The  final  column  of  the  table,  Indeps  Call,  refers  to  a  command  that  is  required 
before  a  table  in  the  file.  This  command  defines  the  independent  variables  to  be  used  to 
define  the  state  of  the  corresponding  table/property.  Every  parameter  marked  N/A  is 
handled  by  the  first  indeps  call  that  specifies  enthalpy  and  pressure  as  the  two 
independent  variables  by  which  the  tables  are  organized.  The  four  parameters  that  have 
the  indeps  call  field  in  the  preceding  table  require  an  additional  call  before  that  particular 
parameter  table  is  instantiated.  Note  the  case  sensitive  nature  of  the  command,  enthalpy 
as  a  function  of  temperature  and  pressure  has  a  lower  case  “i”  in  the  indeps  command. 
The  format  of  the  call  is  described  in  the  NPSS  thermodynamics  reference  sheets  [22] 
and  an  example  is: 

hTindeps  =  {"Tt",  "Pt"}; 

Enthalpy  as  a  function  of  temperature  and  pressure  was  the  easiest  table  to  create 
as  it  is  simply  a  refonnatting  of  the  raw  NIST  data.  Data  was  collected  from  NIST  and 
loaded  into  an  Excel  spreadsheet.  In  the  spreadsheet,  the  data  was  organized  by  the 
independent  variables  and  then  it  was  wrapped  into  ModelCenter  where  the  spreadsheet 
contents  could  be  loaded  into  arrays.  The  arrays  were  then  sent  to  a  Visual  Basic  script 
where  they  were  organized  into  the  NPSS  table  fonnat  and  output  to  a  text  file.  The  file 
was  saved  for  future  integration  into  the  overall  methane  .fpt  file. 

The  remaining  tables  divide  into  two  groups,  enthalpy  based  data  and  entropy 
based  data.  Two  Kriging  estimators  were  built  based  on  those  input  variables.  The 
process  for  building  a  Kriging  model  is  to  save  the  Excel  data  as  tab  delimited  text  data 
for  easy  import  into  ModelCenter  using  the  data  import  function.  Following  the 
ModelCenter  data  import  dialogue  allows  you  to  specify  delimiter  method,  column 
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headings,  columns  for  import,  and  input/output  status  of  a  parameter.  Once  the  data  is 
loaded  with  the  appropriate  specified  input  parameters  (enthalpy  and  pressure  or  entropy 
and  pressure),  the  RSM  toolkit  is  launched  so  a  Kriging  model  can  be  generated.  The 
best  results  for  the  Kriging  come  from  an  even  spacing  of  data  over  the  entire  range  of 
interest  and  fewer  than  1000  data  points,  which  causes  a  tradeoff  between  data  density 
and  range  coverage.  The  RSM  toolkit  provides  some  basic  measures  of  quality  for  the 
Kriging  estimator,  an  adjusted  R"  and  a  1  to  5  star  rating.  The  first  attempt  at  Kriging 
resulted  in  less  than  adequate  results. 

Experimentation  with  the  Kriging  estimator  led  to  a  more  complicated  solution 
that  reduced  error  and  somewhat  reduced  the  impact  of  the  tradeoff  between  data  density 
and  range  coverage.  The  data  was  divided  by  phase.  This  created  four  Kriging 
estimators  to  generate  the  remaining  tables:  the  division  of  liquid  versus  supercritical  and 
the  basis  of  enthalpy  or  entropy  as  the  independent  variable.  Additional  data  was  taken 
for  a  more  complete  coverage  of  the  range,  all  of  the  data  was  reorganized  by  phase,  and 
the  four  Kriging  estimators  were  generated.  The  change  of  the  organization  of  the  data 
from  temperature  to  enthalpy  as  the  independent  variable  caused  some  data  in  the  range 
of  the  vapor  phase  to  be  included  in  the  input  data  (the  higher  enthalpy  input  combined 
with  the  lower  pressure  input).  No  data  in  this  range  was  included  in  the  generation  of 
the  estimators  so  it  is  assumed  that  the  Kriging  provides  incorrect  data.  Therefore, 
additional  logic  was  included  to  eliminate  the  vapor  phase  data  from  the  fluid  property 
tables.  The  logic  is  explained  below  in  the  Rejection  of  Vapor  Data  section  as  the  same 
method  was  included  in  the  oxygen  tables. 
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Enthalpy,  entropy,  and  pressure  input  values  are  then  organized  into  ModelCenter 
arrays  and  fed  into  the  Kriging  estimators  in  loops  to  get  complete  coverage  of  the  input 
data  ranges.  Examining  the  first  combination  of  enthalpy-based  liquid  data  Kriging 
estimator  as  an  example,  the  data  generation  process  consists  of  an  outer  enthalpy  loop 
with  a  nested  pressure  loop.  This  organization  of  data  generation  is  based  around  the 
required  NPSS  table  organization.  In  other  words,  for  each  input  enthalpy,  the  inner  loop 
calculates  the  Kriging  estimate  of  each  thermodynamic  parameter  for  each  specified 
pressure.  The  output  from  the  Kriging  estimator  is  a  single  value  of  each  thermodynamic 
property  of  the  data  type  double.  The  output  doubles  are  stored  in  arrays  of  doubles  by 
the  “Pressure  For  Each  Loop”  component.  Because  ModelCenter  passes  arrays  from 
component  to  component  in  the  form  of  arrays  of  strings,  some  additional  scripting  was 
required  to  organize  the  data  before  preparing  it  for  output.  A  Javascript  component 
converts  the  array  of  doubles  into  a  string.  The  single  string  for  each  property  includes  n 
values  of  each  parameter  where  n  is  the  number  of  input  pressures.  The  enthalpy  For 
Each  Loop  then  collects  each  string  into  an  array  of  strings.  Therefore,  after  all  input 
data  is  processed  through  the  Kriging  estimator,  for  each  parameter,  an  array  of  strings 
exists,  one  string  for  each  value  of  input  enthalpy.  Each  string  contains  the  n  parameter 
values  calculated  based  on  the  input  pressure  values. 

This  process  is  repeated  four  times,  once  for  each  Kriging  estimator.  The  result  is 
four  string  arrays  that  are  already  organized  in  the  manner  required  by  the  NPSS  table 
syntax.  The  data  has  to  be  organized  monotonically  within  a  table  for  all  input  variables. 
Care  must  be  taken  to  avoid  overlapping  enthalpy  values  when  using  multiple  Kriging 
estimators  for  one  table.  The  phase  change  boundary  provided  a  straight  forward 
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demarcation  of  the  input  data,  but  overlap  and  resulting  non-monotonic  tables  could  be  an 
issue  when  creating  multiple  Krigings  without  the  clear  demarcation.  An  output  file  is 
generated  using  the  proper  syntax  and  the  string  data  is  placed  properly  in  the  table. 

The  total  gas  constant  is  a  constant  for  the  pure  chemical  species  and  was 
therefore  provided  as  a  function  rather  than  a  table.  The  previously  developed  enthalpy 
as  a  function  of  pressure  and  temperature  table  must  be  added  to  the  final  fluid  property 
table  file.  A  future  upgrade  could  include  the  integration  of  the  different  table  generation 
techniques. 

Oxygen  FPT  Generation 

The  Oxygen  fluid  property  table  generation  process  was  very  similar  to  the 
methane  table  generation  except  for  a  few  complications.  First,  more  data  was  available, 
so  the  temperature/enthalpy  data  range  is  significantly  larger  than  that  of  the  methane 
table.  Additionally,  the  constant  pressure  specific  heat  as  a  function  of  enthalpy  and 
pressure  in  the  super  critical  phase  Kriging  estimator  was  of  low  quality.  The  ratio  of 
specific  heats  was  predictably  also  poor  (see  equation  38).  Temperature  data  was 
generated  based  on  enthalpy  and  pressure  input  into  the  Kriging  estimator  and  the 
definition  of  constant  pressure  specific  heat  was  employed  via  a  numerical  derivative. 


The  numerical  derivative  was  taken  by  sorting  the  oxygen  NIST  data  by  pressure, 
and  then  calculating  the  partial  derivative  via  equation  37  [24]. 


Cp\x 


II 


x+l 


lx+l 


Hx-i\ 
-  Tx- 1 ) 


(37) 
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The  Cp  values  at  the  edges  of  each  grouping  of  constant  pressure  data  were  then 
queried  and  supplied  directly  from  the  NIST  database  [23].  Data  points  near  the  critical 
point  and  near  the  boiling  point  were  also  directly  supplied  from  NIST  to  reduce  error. 
The  Excel  file  was  then  wrapped  in  ModelCenter  to  save  the  different  columns  of  data 
into  arrays.  The  arrays  were  then  output  to  text  in  the  appropriate  NPSS  table  syntax. 

The  ratio  of  specific  heats  was  simply  included  in  the  fluid  property  table  as  a  function  by 
its  definition: 

y  =  /cv  (38) 

When  gamma  is  queried  by  NPSS,  the  Cp  and  Cv  at  the  given  state  will  be  found 
in  the  table  and  gamma  will  be  calculated  directly.  The  oxygen  gamma  function  is  listed 
below  as  an  example  of  the  format  for  a  NPSS  fluid  property  table  function. 

real  gam  (real  ht,  real  Pt)  {  return  (Cp(ht,Pt)/Cv(ht,Pt)); 

Oxygen  tables  from  all  three  sources  must  then  be  manually  combined  into  one 
fluid  property  table  file.  The  input  sources  are  the  master  table  created  from  the  four 
Kriging  estimators,  the  enthalpy  as  a  function  of  temperature  and  pressure  table  generated 
directly  from  the  NIST  data,  and  the  independently,  spreadsheet  generated  constant 
pressure  specific  heat  table. 

Rejection  of  Vapor  Data 

NIST  data  on  fluid  properties  is  supplied  as  a  function  of  temperature  and 
pressure.  NPSS  needs  the  data  organized  as  a  function  of  enthalpy  and  pressure.  In 
creating  the  fluid  property  tables,  NIST  data  was  captured  and  the  vapor  phase  of 
methane  and  oxygen  were  excluded  prior  to  Kriging  estimator  generation.  After  the  data 
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was  reorganized  by  an  enthalpy  basis,  the  input  ranges  of  enthalpies  and  pressures  were 
an  easily  organized  matrix  of  values  based  on  phase  boundaries,  executed  with  nested  for 
loops.  This  presents  a  problem  because  the  new  enthalpy  pressure  input  data  ranges 
include  regions  that  are  vapor,  and  additionally,  because  of  the  exclusion  of  vapor  data  in 
the  original  sampling,  the  data  is  rendered  meaningless  as  it  does  not  represent  the  reality 
of  the  phase  transition.  The  solution  was  to  eliminate  data  outside  of  the  ranges  of  the 
phase  boundaries.  The  Antoine  equation  allowed  for  the  elimination  of  any  vapor  phase 
data  where  P  is  the  vapor  pressure  in  bar,  T  is  temperature  in  kelvin,  and  A,  B,  and  C  are 
the  Antoine  parameters.  Given  a  temperature,  the  liquid  will  have  a  corresponding  vapor 
pressure. 


log10(P )  =  A  - 


(39) 


The  vapor  pressure  can  also  represent  the  pressure  at  which  a  fluid  of  temperature 
T  will  boil;  therefore,  this  can  be  used  as  a  relationship  between  pressure  and  boiling 
temperature.  The  output  temperature  of  the  Kriging  functions  was  compared  to  the 
Antoine  functions  of  the  respective  fluid.  If  the  temperature  was  higher  than  the  Antoine 
boiling  temperature  for  a  given  pressure,  the  data  was  not  retained  for  tabulation.  In  this 
way,  vapor  phase  data  in  the  new  basis  ranges  was  excluded  from  the  enthalpy  and 
entropy  based  fluid  property  tables. 


Table  10:  Antoine  Coefficients  for  Oxygen  and  Methane. 


Species 

P  -  Unit 

Temp  Range  (K) 

A 

B 

C 

Source 

Oxygen 

Bar 

54.36-  154.33 

3.9523 

340.024 

-4.144 

NIST  [25] 

Methane 

Bar 

90.99-  189.99 

3.9895 

443.028 

-0.49 

NIST  [26] 
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Combustion  Products  FPT  Generation 


The  fluid  property  table  for  the  combustion  products  is  generated  using  NASA’s 
chemical  equilibrium  with  applications.  The  fluid  property  table  was  generated  by 
controlling  CEA  with  ModelCenter  independent  of  NPSS  because  of  the  limitations  of 
the  unreacted  fluid  stations  used  in  the  rocket  model  elements.  Because  the  oxygen  to 
fuel  weight  ratio  is  a  top  level  explored  parameter,  the  combustion  products  fluid 
property  table  requires  three  independent  variables  (pressure,  temperature,  and  OF)  to 
define  the  combustion  products  thermodynamic  state.  Additionally,  a  reverse  lookup 
table  of  temperature  as  a  function  of  enthalpy,  pressure,  and  OF  is  required  per  the  NPSS 
Rockets  Supplement  [18]. 

NPSS  thermochemical  documentation  is  somewhat  vague  about  how  the  program 
would  interpret  the  data  or  whether  O/F  could  even  be  used  by  the  fluid  station  (O/F  is 
not  listed  as  a  property  supported  by  fluid  property  tables  in  the  documentation,  although 
it  is  a  variable  which  can  be  queried  at  a  fluid  station).  Therefore,  it  was  decided  to 
generate  a  unique  combustion  products  fluid  property  table  for  every  iteration  of  O/F. 

The  appropriate  FPT  would  have  been  generated  before  NPSS  calculations  and  then 
called  from  NPSS  by  using  ModelCenter  to  update  the  NPSS  run  file  based  on  the  O/F 
input.  This  solved  the  problem  with  the  documentation  ambiguity  and  should  have 
reduced  error  by  limiting  the  algorithm  to  two  interpolations  (temperature  and  pressure) 
rather  than  three  interpolations.  This  method  was  abandoned  when  an  example  fluid 
property  table  for  hydro  gen/oxygen  combustion  products  was  generated  using  the  NPSS 
Rocets  database. 
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Ultimately,  the  chosen  method  was  to  generate  a  large  three  dimensional  table 
which  NPSS  would  linearly  interpolate  over  three  independent  variables  for  gas 
properties.  A  ModelCenter  program  generated  CEA  input  files  for  the  rocket  problem, 
and  then  captured  the  CEA  output.  The  captured  data  was  then  organized  and  output  into 
the  proper  NPSS  fluid  property  table  format. 

The  ranges  of  pressure  and  temperature  for  the  combustion  products  FPTs  were 
detennined  from  historical  examples  found  in  an  overview  of  rocket  engine  parameters 
by  Oskar  J  Haidn  [27].  From  the  listed  rocket  engines,  the  upper  stage  engines  ranged  in 
chamber  pressure  from  522  to  2147  psia.  The  overall  range  of  liquid  propellant  rockets 
to  include  booster  and  main  stage  engines  is  522  to  3698  psia.  Additionally,  the  Russian 
RD-167  methane/lox  rocket  engine  generates  a  chamber  pressure  of  about  2500  psia  [6]. 
The  limits  of  the  fluid  property  table  for  the  combustion  products  will  therefore  be  1 5 
psia  to  4000  psia.  Based  on  Hall’s  approximate  highest  chamber  temperature,  7000  R 
[8],  the  range  for  the  combustion  products  temperature  independent  variable  will  be  1000 
R  to  8000  R.  The  step  size  for  both  parameters  will  be  small  to  reduce  interpolation 
error.  The  step  size  for  both  temperature  and  pressure  is  the  driver  behind  the  error 
associated  with  the  portion  of  this  FPT  that  is  based  on  temperature,  pressure,  and  O/F. 

The  O/F  ratios  included  in  this  table  are  based  on  CEA  analysis  using  the  Rocket 
problem.  As  expected,  the  optimum  O/F  ratio  for  methane  and  oxygen  is  at  a 
significantly  different  ratio  than  hydrogen  and  oxygen  and  can  be  seen  in  Figure  6.  CEA 
is  used  as  a  first  look  to  determine  a  neighborhood  of  interest  for  the  O/F  parameter.  The 
resulting  vacuum  Isp  is  calculated  based  on  only  chemical  and  thermodynamic  properties 
and  has  no  system  level  dependence  on  energy  or  material  transport.  It  is  therefore  only 
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used  to  examine  relationships.  Additionally,  the  calculated  Isp  seen  in  Figure  6  is  based 
on  frozen  flow  conditions,  and  therefore  represents  a  low  estimate.  The  low  number  is 
not  a  concern  as  the  input  to  CEA  is  not  based  on  any  optimization. 


Figure  6:  CEA  Isp  vs.  O/F  at  Various  Expansion  Ratios  for  CH4/02,  Frozen  Flow. 


Changing  the  expansion  ratio  has  a  significant  effect  on  the  specific  impulse  and 
the  location  of  the  optimum  O/F  ratio.  Note  the  Hall  hydrogen  concept  had  an  expansion 
ratio  of  4.37.  This  first  look  shows  that  the  area  of  interest  for  the  current  work  should 
include  O/F  ratios  from  2  to  5. 

Temperature  as  a  Function  of  Enthalpy,  Pressure,  and  OF 
Per  the  NPSS  User’s  Guide  [19]  and  the  RocketCombl  Element  [18],  the  fluid 
property  table  requires  a  reverse  lookup  table  to  find  temperature  as  a  function  of 
enthalpy,  pressure,  and  O/F.  The  first  attempt  at  table  generation  used  the  method  that 
was  employed  in  the  propellant  tables,  a  Kriging  interpolation  that  allowed  for  estimates 
of  fluid  properties  to  be  organized  by  the  different  variables.  For  reasons  discussed  in  the 
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results  section,  this  was  unsuccessful  and  error  was  significant.  A  new  approach  was 
attempted  for  this  table,  a  three  dimensional  cubic  polynomial.  Significantly  more  data 
was  used  in  the  generation  of  the  polynomial  as  compared  to  a  Kriging  generation  which 
is  limited  to  1000  data  points.  Additionally,  data  density  was  increased  in  the  range  of 
interest  where  convergence  is  expected.  The  final  polynomial  employed  in  the  model  is 
in  Table  11. 


Table  11:  Polynomial  Approximation  of  Temperature  as  a  Function  of  H,  P,  and  OF. 


Terms 

Coefficients 

Constant 

6.6934E+02 

O/F 

3.4226E+03 

Enthalpy 

1.5122E+00 

Pressure515  Pressure 

-3.068  IE-04 

O/F*  O/F 

-6.6019E+02 

Enthalpy*Enthalpy 

-1.2292E-04 

Pressure*0/F 

5.2064E-01 

Pressure*Enthalpy 

2.7619E-04 

0/F*Enthalpy 

-5.1952E-01 

Pressure  *  Pressure  *  Pressure 

4.8916E-08 

O/F*  O/F*  O/F 

3.5103E+01 

Enthalpy*Enthalpy*Enthalpy 

4.9080E-09 

Pressure*Pressure*0/F 

-2.5999E-05 

Pressure  *  Pressure  *  Enthalpy 

-3.8899E-08 

Pressure  *  O/F  *  O/F 

-5.51 14E-02 

Pressure  *  O/F  *  Enthalpy 

-1.7437E-05 

Pressure*Enthalpy*Enthalpy 

-1.1080E-08 

O/F  *  O/F  *  Enthalpy 

6.1895E-02 

0/F*Enthalpy*Enthalpy 

9.2284E-06 

This  function  is  simply  placed  in  the  fluid  property  table  in  the  proper  fonnat  and  returns 
a  temperature  when  queried.  The  format  is  similar  to  the  more  simplistic  gamma 
function  described  above. 


57 


Thermodynamic  Reference  State 

Each  source  of  data  is  defined  by  a  different  thermodynamic  reference  state.  The 
NIST  data  defines  the  thennodynamic  reference  state  as  internal  energy  is  equal  to  zero  at 
273.16  K  for  the  saturated  liquid  [25,  26].  This  also  represents  the  enthalpy 
thennodynamic  reference  state  per  the  definition  of  enthalpy  [28]: 

H  =  U  +  pV  (40) 

This  is  true  for  each  propellant.  The  combustion  products  data  set  is  based  on  the 
CEA  thennodynamic  reference  state,  which  is  not  directly  and  simply  defined.  For  the 
purposes  of  the  MDEAN,  this  is  not  necessary.  The  reference  value  is  arbitrary  because 
the  change  in  enthalpy  from  one  thermodynamic  state  to  another  is  the  important  relation, 
and  by  the  definitions  of  the  thennochemical  properties,  a  change  from  one 
thermodynamic  state  to  another  will  represent  the  same  change  in  enthalpy,  internal 
energy,  or  entropy  regardless  of  the  starting  value.  The  difference  in  thermodynamic 
reference  state  amongst  the  different  data  sets  is  unimportant  except  when  the  chemicals 
of  the  different  states  interact.  The  only  interaction  takes  place  in  the  combustion 
chamber.  The  NPSS  RocketCombl  element  already  takes  the  reference  state  into 
account  via  the  correction  term  Href  used  when  calculating  the  htMix,  the  enthalpy  of 
mixture  of  propellants  as  they  enter  the  combustion  chamber. 

Hmix  =  (Hf  -  Hf  ref)Wf  +  ( Hox  -  Hoxref)Wox  (41) 

The  values  for  the  thermodynamic  reference  state  correction  terms  are  calculated 
based  on  the  available  CEA  data.  The  thennochemical  library  included  with  CEA  has 
one  data  point  for  each  oxygen  and  methane  [20].  The  NIST  database  was  then  queried 
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at  the  same  pressure  and  temperature  and  the  density  and  enthalpy  were  recorded.  The 
data  can  be  seen  in  Table  12. 


Table  12:  CEA  and  NIST  Propellant  Thermodynamic  Reference  Data. 


CEA 

BP  (K) 

P  (atm) 

Density  (g/mL) 

Enthalpy  (KJ/mol) 

CH4 

111.643 

1 

0.4211 

-89.233 

02 

90.17 

1 

1.149 

-12.979 

NIST 

CH4 

111.64 

1 

0.42239 

-0.0013586 

02 

90.17 

1 

1.1413 

-4.2686 

Each  Href  can  then  be  calculated  for  each  propellant  by  taking  the  difference 
between  the  two  respective  enthalpies  and  then  converting  to  the  NPSS  unit  of  enthalpy, 
Btu/lbm. 


K]  Btu 

Hf  ref  =  (-0.0013586  - 89.233)  =  89.2316414— ^  =  2391.6893 

1,1  mol 


lbr 


(42) 


KJ  Btu 

Hnxref  =  (-4.2686  - 12.979)  =  8.7104—! -  =  117.0293 

'  J  m  n  I 


mol  lbm  (43) 

The  values  of  the  reference  enthalpy  conversions  are  entered  into  the  respective 
fuel  streams  in  NPSS  in  the  section  for  entering  initial  guesses. 


DEAN  Changes  to  Accommodate  Methane  Fuel 

Several  updates  and  changes  were  made  to  the  DEAN  model  to  achieve  a 
methane  based  simulation.  The  starting  point  of  the  MDEAN  was  the  sixth  iteration  of 
the  hydrogen  DEAN  developed  by  Simmons.  Simmons’  work  was  based  on  Hall’s  work, 
and  Hall’s  work  was  built  on  the  original  DEAN  NPSS  produced  by  Martin. 

Development  of  the  sixth  and  seventh  versions  of  the  hydrogen  based  DEAN  continued 
in  parallel  with  the  development  of  the  MDEAN. 
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Conical  Approximation 

As  stated  in  the  background  section,  a  large  departure  from  Hall’s  work  is  the 
reliance  on  conical  approximations  with  regards  to  the  aerospike  contour.  Hall’s  DEAN 
work  demonstrated  that  the  conical  aerospike  developed  reasonable  approximations  for 
mass  calculations  while  significantly  decreasing  processing  time  [8].  Therefore,  all 
nonlinear  code  was  removed  from  the  DEAN  code  in  the  transition  to  the  MDEAN. 

Pre-Processing 

Work  by  Simmons  on  the  sixth  version  of  DEAN  updated  the  linkage  of  top  level 
design  parameters  to  underlying  model  behavior.  A  major  limitation  with  the  existing 
DEAN  was  the  length  of  time  required  for  examining  the  input  variable  trade  space.  It 
was  possible  and  very  likely  to  explore  a  rocket  concept  that  was  not  physically 
realizable.  DEAN6  attempts  to  link  internal  NPSS  geometry  inputs  into  well  understood 
top  level  geometric  parameters,  thereby  greatly  reducing  the  number  of  physically 
meaningless  design  points.  At  the  same  time,  this  provides  the  DEAN  a  clear  set  of  input 
parameters  that  can  be  explored  via  the  ModelCenter  parametric  study  function  and  the 
design  of  experiments  function.  The  input  variables  are  listed  in  Table  13. 

Table  13:  Methane  DEAN  Design  Variables. 


Design  Variable 

Unit 

Expansion  Ratio 

N/A 

Throat  Area 

in2 

Chamber  Length 

in 

Characteristic  Length 

in 

Thrust 

lbf 

O/F  Ratio 

N/A 

From  these  input  variables,  several  important  geometric  values  which  define  a 
specific  rocket  design  iteration  can  be  calculated.  First,  used  as  direct  input  to  the  NPSS 
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rocket,  are  the  chamber  and  throat  outer  and  inner  radii.  These  four  values,  combined 
with  Martin’s  original  NPSS  assumptions  about  the  shape  of  an  aerospike  engine  define 
the  contour  that  combines  the  combustion  chamber  to  the  converging  diverging  nozzle. 
The  original  geometry  of  the  chamber  is  split  up  into  6  stations,  5  at  the  midpoint  of  equal 
lengths  of  the  chamber  and  the  final  station  at  the  throat.  Station  five  and  six  are  in  the 
converging  portion  of  the  design  where  the  outer  chamber  wall  slopes  inward.  At  station 
six,  the  chamber  wall  radius  is  one  inch  narrower  than  the  chamber  radius  at  the  injection 
face  and  station  five  is  half  way  up  the  linear  contraction  at  0.5  inches  narrower  than  the 
chamber  radius  at  the  injection  face.  A  ModelCenter  rendered  aerospike  contour  is 
included  in  Figure  7  for  clarity.  Note  that  this  is  the  initial  estimate  of  conical  geometry 
used  to  derive  relationships  required  by  NPSS.  Final  geometry  uses  a  more  realistic 
approximation  for  volume  and  mass  calculations. 


Figure  7:  DEAN  Conical  Aerospike  Contour. 


In  addition  to  the  four  radii  that  define  the  aerospike,  the  initial  geometric 

•  •  .  .  ^ 
calculations  detennine  the  volume  of  the  chamber  based  on  the  characteristic  length,  L  , 
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which  is  equal  to  the  chamber  volume  divided  by  the  area  of  the  throat.  The  volume  of 
the  chamber  is  also  an  input  required  by  NPSS.  For  a  quick  check  on  the  feasibility  of 
the  input  geometry,  the  injection  ratio  is  calculated  as  the  ratio  of  the  chamber  area  at  the 
injector  face  to  the  area  of  the  throat.  For  the  input  geometry  to  be  reasonable,  this  value 
must  be  greater  than  about  two  [4].  This  limit  is  based  on  the  areas  of  the  throats  that  are 
likely  to  converge  given  the  requirement  of  25,000  lbf  of  thrust.  This  is  an  imperfect  tool 
as  the  literature  relates  the  injection  ratio  to  the  diameter  of  a  circular  throat,  but  the 
constraint  provides  a  useful  rule  of  thumb  that  eliminates  the  possibility  of  losing  the 
converging  portion  of  the  converging/diverging  nozzle.  The  final  geometric  calculation 
to  be  input  into  NPSS  is  the  length  of  the  aerospike  nozzle.  The  aerospike  nozzle  length 
is  approximately  the  length  of  a  cone  nozzle  with  a  half  angle  of  12  degrees  [7]. 

Initial  Input  to  Cooling  Channel  Pressure  Profiles 

When  the  MDEAN  fails  to  converge,  it  is  usually  related  to  the  independent 
variables  and  the  dependent  relationships  in  the  propellant  cooling  channels.  In  an 
attempt  to  increase  the  robustness  of  the  model,  the  initial  guesses  for  the  pressure  profile 
through  the  cooling  channels  were  tied  to  input  parameters.  Previous  iterations  of  the 
DEAN  relied  on  hard  numbers  coded  into  the  input  guesses,  presumably  captured  from 
converged  cases.  The  method  of  developing  the  oxygen  pressure  profile  begins  with  the 
input  design  variables  throat  area  (At)  and  thrust  (F).  Additionally,  a  guessed  value  of  Isp 
is  required.  For  the  MDEAN,  this  guess  is  350  seconds.  These  input  values  and  the 
rough  guess  calculate  the  mass  flow  rate: 


With  a  guess  for  c*,  the  characteristic  exhaust  velocity,  or  6020  ft/s  [7],  a  rough 
estimate  for  chamber  pressure,  Pc,  can  be  calculated  with  the  mass  flow  rate  and  the  area 
of  the  throat,  At: 

me* 

Pc=^~  (45) 

Now  the  pressure  is  known  at  either  ends  of  the  propellant  feed  system.  An 
approximately  5%  drop  in  pressure  is  assumed  from  the  tanks  to  the  pump  inlets.  The 
pressure  ratios  for  all  three  pumps  are  input  parameters,  so  the  pressure  profiles  up  to  the 
cooling  jacket  can  be  estimated.  Starting  at  the  other  end  of  the  feed  system,  a  20% 
pressure  loss  is  assumed  for  the  injector.  The  pressure  drop  across  the  turbine  is 
estimated  by  taking  an  average  over  several  converged  cases  and  resulted  in  a  25%  loss 
for  the  fuel  and  a  50%  loss  for  the  oxidizer.  The  oxidizer  turbine  bypass  line  was 
calculated  as  the  average  of  the  entrance  and  exit  of  the  oxidizer  turbine.  Therefore,  for 
both  propellant  feed  systems,  the  pressure  profile  has  been  estimated  up  to  the  entrance 
and  past  the  exit  of  the  cooling  channels,  providing  a  change  in  pressure  over  the  cooling 
channels.  This  pressure  drop  is  modeled  as  a  constant  pressure  drop  across  the  stations. 
These  calculations  provide  an  estimate  of  the  pressure  profile  that  is  flexibly  based  on  the 
input  parameters  and  therefore  increases  the  robustness  of  the  model. 

NPSS  Updates 

The  specific  changes  to  the  NPSS  code  are  outlined  in  Appendix  B.  The  goals 
and  implications  of  the  updates  are  discussed  in  this  section. 
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Thermodynamic  Updates 

The  first  major  change  is  the  command  to  utilize  the  new  sources  of 
thermodynamic  properties.  At  the  beginning  of  any  NPSS  model,  if  the  standard 
databases  will  not  be  used,  an  alternate  must  be  set  with  the  SetThermoPackage 
command.  As  discussed  previously,  the  MDEAN  uses  the  fluid  property  table 
thermodynamic  package.  To  use  the  developed  tables,  the  composition  of  each  flow  has 
to  be  set  to  match  the  fluid  property  table  file  name  and  this  is  usually  set  as  the  comp 
variable  within  the  instantiation  of  an  element.  By  using  the  SetDefaultComposition 
command  outside  an  element,  any  non-specified  flow  station  compositions  are  set  to  the 
specified,  default  value.  The  organization  of  the  DEAN  NPSS  model  makes  it  easy  to  set 
the  proper  compositions.  In  the  DEAN  model,  methane  must  be  set  as  the  default 
composition  for  the  methane  feed  system  and  must  be  specified  as  the  fuel  inlet 
composition  in  the  RocketCombl  (the  combustion  chamber)  element.  Similarly,  oxygen 
is  set  at  the  beginning  of  the  oxygen  feed  system  as  the  default  and  is  specified  at  the 
combustion  chamber  element  oxidizer  inlet.  The  combustion  products  are  set  as  the 
default  composition  before  the  combustion  chamber  and  are  specifically  set  at  the 
combustion  chamber  exit. 

As  discussed  earlier,  the  thermodynamic  bases  must  match  amongst  the  different 
species/mixtures.  This  is  accomplished  through  calculations  in  the  combustion  chamber 
element,  but  the  thennodynamic  reference  states  must  be  specified  in  the  body  of  the 
MDEAN  code.  Both  COMB.Ful.htRef  and  COMB.Fl  oxid.htRef  are  specified  in  the 
MDEAN  section  of  code  that  outlines  initial  values  to  start  the  solver  iteration. 
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Cooling  Channel  Enthalpy  Profile 

Both  the  fuel  and  oxidizer  have  guesses  for  the  enthalpy  profiles  through  the 
respective  feed  systems.  The  original  values  are  assumed  to  be  based  on  converged 
solutions.  The  profile  for  each  species  was  converted  to  the  NIST  thennodynamic  basis. 
This  was  accomplished  by  starting  the  calculation  update  at  the  pumps,  where  the 
thennodynamic  state  is  defined  by  pressure  and  temperature.  The  existing  NPSS 
assumed  increases  in  enthalpy  were  then  applied  through  the  flow  system.  Several 
failures  of  convergence  reference  the  enthalpy  of  the  different  cooling  channels.  Each 
cooling  channel  enthalpy  is  an  independent  variable  which  is  perturbed  in  the  solver.  The 
methane  side  enthalpy  profile  was  later  replaced  by  data  from  a  converged  run.  Future 
work  in  the  DEAN  or  MDEAN  should  include  the  upgrade  of  this  system  to  tie  these 
initial  enthalpy  guesses  to  model  input  values.  The  updated  values  for  pressure  and 
enthalpy  lead  to  a  new  defined  thermodynamic  state  and  therefore  the  density  estimates 
for  the  same  flow  stations  were  updated. 

Pump  Changes 

In  previous  DEAN  work,  the  pressure  ratios  of  the  fuel  pumps  were  set  equal  to 
each  other  as  targets  for  convergence.  This  was  done  by  creating  a  dependent  condition 
in  the  top  level  DEAN  of  the  pressure  ratio  of  the  first  pump  is  equal  to  the  pressure  ratio 
of  the  second  pump.  This  usually  resulted  in  a  small  pressure  rise  from  the  first  pump, 
followed  by  a  large  pressure  rise  in  the  second  pump.  This  was  updated  in  the  sixth 
version  of  the  DEAN  so  that  each  fuel  pump  produced  an  equal  pressure  rise.  This  was 
accomplished  by  adding  the  variable  deltaP  to  the  pump  element.  The  dependent 
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condition  in  the  DEAN  was  then  updated  so  the  variable  deltaP  was  converged  over  the 
two  fuel  pumps.  The  DEAN  now  requires  the  custom  element,  Pump02.int. 

The  fuel  pump  rotational  rate  in  RPM  required  update  based  on  the  fluid 
properties  of  the  new  fuel.  The  goal  for  the  update  was  to  maintain  the  specific  speed, 
Ns,  of  both  fuel  pumps.  The  relation  uses  U.S.  customary  units  as  specified  [7]: 


where 

N  =  Actual  Pump  Speed  (RPM) 

Qe  =  Volumetric  flow  (ft3/s) 

He  =  Pump  Head  (ft) 

By  maintaining  specific  speed,  the  same  class  of  pump  is  preserved;  therefore  the 
estimates  used  for  the  pump  efficiencies  can  be  maintained.  The  hydrogen  pumps 
operated  at  a  speed  of  1 10,000  RPM  and  this  was  updated  to  50,000  RPM  for  the  more 
dense  methane  based  on  a  different  fuel  flow  rate  (the  O/F  is  significantly  different),  the 
change  in  the  nature  of  the  fuel  pump  pressure  ratios,  and  the  different  fuel  density. 

Constraints 

The  DEAN  model  can  converge  onto  designs  that  are  physically  unrealizable. 

The  constraints  section  of  the  DEAN  attempts  to  automatically  analyze  and  eliminate 
converged  models  that  violate  a  fundamental  concept  of  operation.  Currently,  the 
constraints  section  of  code  checks  for  supersonic  flow  of  propellant  in  the  cooling 
channels  and  the  combustion  chamber.  Supersonic  flow  in  these  regions  is  inappropriate 
as  it  physically  unrealizable  without  a  converging  diverging  nozzle.  Sonic  flow  in  these 
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regions  can  occur  through  an  expansion  of  boundary  layers.  This  represents  an 
undesirable  design  because  it  involves  an  inefficient  high  pressure  drop  through  the  flow 
system  and  can  cause  shocks  in  the  flow. 

Cooling  Channel  Mach  Number 

An  important  constraint  in  the  development  of  a  mathematically  reasonable  model 
that  is  also  physically  realizable  is  the  elimination  of  supersonic  fluid  velocities  anywhere 
in  the  design  other  than  the  rocket  nozzle.  The  propellant  cooling  channels  constrict  to 
small  cross  sectional  areas  to  increase  heat  transfer  and  accommodate  the  chamber  and 
nozzle  geometry.  By  simple  incompressible  flow,  the  smaller  flow  area  will  accelerate 
the  fluid.  Additionally,  heat  is  flowing  into  the  fluid  from  the  combustion  chamber, 
increasing  the  overall  energy  of  the  propellant.  The  highest  fluid  velocities  will  therefore 
be  found  in  the  cooling  jacket  as  energy  will  dissipate  from  the  fluid  after  the  cooling 
jacket  as  the  two  fluids  are  expanded  through  the  two  turbines. 

The  Mach  numbers  must  stay  within  limits  that  would  cause  shocks  in  the  flow. 
For  oxygen,  the  Mach  number  must  stay  below  0.6.  For  methane,  the  Mach  number  must 
remain  below  0.9  [13].  The  sonic  velocity  for  the  propellants  is  a  function  of  fluid  state 
as  defined  by  temperature  and  pressure  in  this  case.  Again,  NIST  data  was  available  over 
the  range  of  interest  for  both  fluids.  To  automate  the  process  of  checking  whether  a  run 
remained  within  the  required  boundary,  a  Kriging  interpolation  was  used  to  estimate  the 
sonic  velocity  of  the  fluid  given  the  temperature  and  pressure  of  the  fluid  station  as 
reported  by  NPSS.  After  NPSS  has  converged,  for  every  fluid  station  in  the  cooling 
jacket,  the  Mach  number  is  calculated  by  dividing  the  fluid  velocity  by  the  Kriging 
estimate  of  the  sonic  velocity  and  if  the  largest  Mach  number  is  greater  than  the  above 
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limits,  an  error  occurs,  warning  the  user  of  the  violated  constraint.  Below  is  a  discussion 
of  the  validity  of  the  Kriging  interpolation. 

The  oxygen  Kriging  interpolation  has  an  average  error  of  2.62%  over  201 
observations  that  span  the  temperature  and  pressure  range  of  the  function.  The  Kriging 
design  range  was  developed  by  observations  of  NPSS  behavior  while  remembering  the 
ranges  recorded  by  the  converged  design  points  developed  in  past  DEAN  research. 
Original  estimates  were  insufficient  in  range  as  some  iterations  of  NPSS  use  very  large 
increases  in  pressure.  The  oxygen  Kriging  was  developed  using  a  temperature  range  of 
100  to  1500  Rankin  and  a  pressure  range  of  15  to  5000  psia.  Note  the  Kriging  was 
developed  without  any  vapor  phase  data  and  therefore  if  the  input  state  of  a  fluid  station 
is  vapor,  the  Kriging  will  not  deliver  an  accurate  sonic  velocity.  This  is  acceptable 
because  the  design  of  the  system  assumes  only  liquid  and  supercritical  phases  are  present 
in  the  propellant  feed  systems.  A  comparison  of  NIST  data  to  the  Kriging  interpolation 
of  the  data  at  the  same  state  as  the  given  NIST  data  shows  the  Kriging  data  is  a  good  fit 
and  is  suitable  for  the  automated  process  (see  Figure  8).  The  data  used  for  the 
comparison  of  NIST  versus  Kriging  data  from  which  the  average  error  was  derived  is  a 
different  set  of  NIST  data  from  which  the  Kriging  interpolation  function  was  derived. 

Per  the  Kriging  algorithm,  the  error  for  the  points  from  which  the  Kriging  interpolation  is 
derived  should  be  zero  or  near  zero  [17].  If  a  selected  run  point  has  a  maximum  Mach 
number  close  to  the  limit  of  0.6,  the  oxygen  fluid  velocity  should  be  examined  against  the 
exact  NIST  data  to  verify  the  validity  of  the  geometric  design. 
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Figure  8:  NIST  vs.  Kriging  Sonic  Velocity  for  Oxygen. 

Blue  =  NIST,  Red  =  Kriging  Estimate;  Sonic  Vel  =  f(T,  P) 

The  fluid  velocities  at  the  methane  stations  are  checked  in  the  same  manner  as  the 
oxygen  stations.  The  methane  sonic  velocity  Kriging  interpolation  method  has  an 
average  error  of  0.37%  over  175  observations  across  the  range  of  the  Kriging.  The 
Kriging  was  developed  over  a  temperature  range  of  193  to  1113  Rankin  and  a  pressure 
range  of  100  to  6000  psia.  NIST  data  at  high  temperatures  was  limited  for  methane 
compared  to  oxygen,  but  the  available  range  should  be  sufficient  for  the  methane  feed 
system.  Again,  any  vapor  data  fed  into  the  Kriging  interpolation  will  provide 
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meaningless  results.  A  comparison  of  NIST  source  data  to  Kriging  generated  data  can  be 
found  in  Figure  9. 


Figure  9:  NIST  vs.  Kriging  Sonic  Velocity  for  Methane. 

Blue  =  NIST,  Red  =  Kriging  Estimate;  Sonic  Vel  =  f(T,  P) 

Chamber  Mach  Number 

Developed  by  Simmons  for  DEAN6,  the  method  for  checking  the  Mach  number 
in  the  chamber  is  based  on  an  analysis  of  the  contraction  ratio,  defined  [4]: 
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Ac  =  Chamber’s  Average  Cross-Sectional  Area  (in") 

At  =  Area  of  the  Throat  (in') 

M  =  Mach  number  in  the  Combustion  Chamber 
y  =  Isentropic  Parameter 

The  logic  in  the  check  subtracts  the  left  side  of  the  equation  from  the  right  hand 
side  to  find  an  error  term  and  then  iterates  the  Mach  number  through  a  ModelCenter 
optimization  tool  until  the  error  term  is  driven  to  zero.  The  resulting  combustion 
chamber  Mach  number  is  displayed  and  should  be  low  to  provide  time  for  combustion. 
The  literature  suggests  the  chamber  Mach  number  should  be  between  0.2  and  0.4  [4], 

Post  Processing 

The  post  processing  calculations  were  largely  inherited  from  the  work  of  Hall. 
Additional  work  to  integrate  the  post  processing  calculations  was  undertaken  by 
Simmons  in  the  development  of  DEAN6. 

Angelino  Contour 

The  MDEAN  NPSS  code  calculates  energy  balances  over  the  system  under  the 
assumption  of  a  bell  nozzle  with  an  expansion  ratio  given  as  an  input.  This  bell  nozzle 
expansion  ratio  can  be  equated  to  an  equivalent  aerospike  nozzle.  Hall  employed  the 
Angelino  method  to  detennine  a  highly  accurate  exit  Mach  number  and  then  a 
corresponding  aerospike  nozzle  length.  This  length  is  then  used  to  develop  the  spike 
non-linear  contour  which  unfortunately  breaks  down  for  calculating  the  geometry  near 
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the  throat.  At  this  point,  Hall  used  the  Angelino  length  in  another  program,  TDK04,  to 
develop  a  working  contour.  Simmons’  work  in  preprocessing  first  estimates  an  aerospike 
length  using  relations  found  in  the  literature  [7],  a  cone  using  a  half  angle  of  12  degrees. 
This  matches  the  Angelino  approximation  well,  but  for  continuity,  the  Angelino 
approximation  is  used  for  the  contour  in  DEAN6  to  avoid  the  use  of  TDK04,  which  is 
now  a  superfluous  step.  Simmons  therefore  adjusted  the  Angelino  curve  to  correct  for 
the  behavior  at  the  throat.  The  Angelino  contour  is  then  incorporated  into  the  next 
section  of  code  which  calculates  the  final  geometry  and  mass  of  the  engine.  When 
comparing  the  Angelino  approximation  to  Hall’s  work,  based  on  TDK04,  the  new 
contour  method  produced  a  significantly  lower  volume  and  therefore  mass.  In  an  effort 
to  maintain  a  conservative  approximation  for  the  thrust  to  weight  ratio,  the  Angelino 
contour  was  averaged  with  the  conical  contour  that  was  calculated  in  the  pre-processing 
steps.  The  resulting  contour  matches  Hall’s  more  sophisticated  contour  calculations  more 
closely  than  the  Angelino  approximation  or  the  conical  approximation.  The  three 
different  contours  can  be  seen  in  Figure  10. 


Conical  Angelino  Averaged 

Figure  10:  Three  DEAN  Estimated  Nozzle  Contours. 
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Engine  Geometry 

Hall’s  first  calculation  detennines  the  combustion  chamber  and  nozzle  internal 
geometry  along  with  the  maximum  temperature  of  the  propellants.  These  are  largely 
defined  through  the  NPSS  code,  but  the  stations  need  update  based  on  the  more  detailed 
aerospike  contour.  The  code  assigns  the  two  nozzle  stations  axial  length  and  radii  to  the 
newly  calculated  non-linear  aerospike  contour. 

These  values  are  then  used  with  material  properties  of  several  candidate  materials 
to  perform  a  structural  analysis.  The  structural  analysis  calculates  wall  thickness  based 
on  the  pressures  found  in  the  NPSS  output.  Additionally,  material  properties  are  supplied 
to  the  ModelCenter  calculations  as  functions  of  temperature.  The  software  applies 
curved  beam  theory  to  detennine  if  the  wall  between  the  propellant  channels  and  the 
combustion  gasses  will  fail.  Based  on  the  NPSS  output  pressures,  the  code  calculates  the 
bending  stress  and  compares  it  to  the  user  specified  method  of  failure:  yield  or  ultimate 
for  the  designated  material.  The  shear  stress  is  calculated  and  must  be  less  than  the 
specified  material  shear  strength,  which  is  defined  here  as  one  third  of  the  ultimate 
strength.  If  failure  occurs,  the  algorithm  increases  the  wall  thickness  iteratively  until  the 
wall  withstands  rocket  operation.  The  wall  thickness  data  and  other  previously  calculated 
geometry  data  are  then  sent  to  the  next  script  for  mass  estimation. 

Thrust  to  Weight  Calculation 

The  final  DEAN  performance  calculation  was  developed  by  Hall  and  integrated 
into  DEAN6  by  Simmons.  The  script  calculates  the  volumes  and  masses  of  the  chamber, 
aerospike,  plumbing  and  the  mass  of  the  turbo-machinery.  The  chamber  and  aerospike 
are  broken  up  into  structural  components  and  cooling  channel  components.  The 
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geometry  is  defined  in  the  previous  ModelCenter  scripts  to  feed  the  volume  calculations 
in  this  script.  With  component  volumes  and  material  densities,  the  mass  of  the  different 
components  are  calculated.  Additionally,  the  turbopump  masses  are  calculated  from  the 
following  engineering  correlation  [4], 

mtp  =  Atb  (48) 

where 

mtp  =  mass  of  the  turbopump  (kg) 

A  =  empirical  coefficient  (1.5  for  development) 
x  =  pump  shaft  torque  (Nm) 

B  =  empirical  exponent  (0.6  for  development) 
Plumbing  length  and  therefore  mass  is  estimated  at  twice  the  engine  length 
(nozzle  length  plus  chamber  length)  while  the  hardware  mass  is  estimated  as  five  percent 
of  the  total  engine  mass.  The  percent  hardware  is  treated  as  a  design  variable  and  can  be 
adjusted  as  necessary.  Five  percent  was  chosen  for  consistency  over  the  different  DEAN 
projects.  All  of  the  mass  is  then  summed  to  find  a  total  engine  mass.  With  the  NPSS 
calculated  thrust,  the  figure  of  performance,  thrust  to  weight,  can  be  calculated. 

The  key  task  in  the  exploration  of  the  MDEAN  is  the  development  of  fluid 
property  tables  that  accurately  provide  thermodynamic  values  for  a  defined  state.  The 
propellant  tables  were  built  using  online  NIST  data  and  the  combustion  products  tables 
were  built  using  the  chemical  property  estimation  program  CEA.  The  tables  were  then 
integrated  to  the  concurrently  developed  DEAN6  model  through  several  changes  to  the 
different  model  components.  After  some  debugging  and  validation,  the  MDEAN  model 
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was  exercised  over  a  wide  range  of  the  input  design  variables.  The  results  of  these  runs 
are  discussed  in  the  following  chapter. 
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IV.  Analysis  and  Results 


This  chapter  will  discuss  the  results  of  the  MDEAN  research.  First,  the  error 
associated  with  the  new  fluid  property  tables  is  explored.  Then  the  results  of  the  runs  of 
the  MDEAN  model  are  examined  to  find  interesting  performance  estimates  over  a  wide 
range  of  the  design  variables.  A  variety  of  thrust  levels  were  explored,  and  one  of  the 
more  promising  design  points  at  25,000  lbf  of  thrust  is  selected  for  a  more  in  depth 
analysis  to  detennine  if  the  rocket  is  physically  realizable.  Finally,  the  chosen  rocket 
concept  is  compared  to  several  other  rockets. 

Oxygen  FPT  Error  Analysis 

The  oxygen  fluid  property  tables  were  more  complex  in  construction  than  the 
methane  fluid  property  tables,  but  were  similar  in  concept  in  that  the  source  NIST  data 
must  be  reorganized  from  being  defined  by  the  independent  variables  of  pressure  and 
temperature  into  data  tables  that  define  the  state  of  the  fluid  by  enthalpy  and  pressure. 
Several  different  methods  of  data  manipulation  were  required  to  obtain  sufficient 
accuracy  in  property  values.  First,  a  Kriging  fit  was  initially  used  over  the  entire  range  of 
required  data  to  include  liquids  and  supercritical  fluids.  Error  for  this  interpolation  was 
such  that  the  resulting  tables  were  deemed  insufficient.  The  data  sets  input  to  the  Kriging 
interpolation  algorithm  generation  were  therefore  split  by  phase.  The  data  boundaries  for 
the  two  Kriging  interpolators  are  found  in  Table  14. 
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Table  14.  Independent  Variable  Ranges  for  Oxygen  Fluid  Property  Tables. 


Phase 

Plow  (psia) 

Pffiah  (psia) 

Hlow  (Btu/lbm) 

HHi2h  (Btu/Ibm) 

Liquid 

15 

4000 

-74 

-11 

Super  Critical 

750 

4000 

-5 

220 

The  quality  of  the  interpolation  algorithm  is  determined  by  simple  comparison  of 
Rf  values  of  the  models,  to  determine  a  qualitative,  but  relative  comparison  of  the  quality 
of  the  estimation.  The  R“  for  the  Cp  Kriging  function  over  the  entire  data  set  was 
72.44%,  insufficient  for  the  purpose  of  tabulating  values.  By  splitting  the  data  sets  by 
phase,  the  liquid  phase  R  was  99.925%,  but  the  super  critical  R  was  still  too  low  for 
use.  To  increase  accuracy  of  the  oxygen  Cp  and  ratio  of  specific  heats  (y) 
thermodynamic  properties  in  the  super  critical  phase  (which  will  be  encountered  in  the 
cooling  channels),  the  following  definitions  of  the  properties  were  used. 


These  definitions  were  employed  in  two  different  ways.  The  y  parameter  is 
simply  solved  for  in  the  NPSS  interpolation  tables  by  replacing  the  table  with  a  function 
that  calls  for  the  corresponding  Cp  and  Cv  values,  which  are  tabulated.  Cp  is  estimated 
through  a  numerical  approximation  of  the  partial  derivative  of  enthalpy  with  respect  to 
temperature  at  constant  pressure.  Temperature,  pressure,  and  enthalpy  data  was 
generated  via  the  super  critical  Kriging  interpolator.  A  sufficient  range  of  enthalpy  and 
pressure  data  was  input,  and  the  temperature  was  estimated  as  an  output.  This  data  was 
placed  in  a  spreadsheet  and  sorted  by  pressure  and  then  the  following  numerical 
derivative  was  applied  to  estimate  the  constant  pressure  specific  heat. 
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(51) 


=  /Hx+i _Jx-A 

V  'x  \  _  rr  I 

WX+1  lX-lJ 

At  the  edges  of  the  blocks  of  constant  pressure  data,  Cp  values  were  taken  directly 
from  the  NIST  database.  Additionally,  data  was  spot  checked  against  the  NIST  database 
and  it  was  found  that  data  near  the  critical  point  (T=278  R,  P=73 1  psia)  had  significant 
error.  Data  in  this  region  was  replaced  with  NIST  data. 

The  fluid  property  tables  were  then  tested  for  accuracy  against  the  NIST  database. 
The  Kriging  estimator  directly  supplies  the  values  that  were  used  in  its  generation  when 
queried  for  an  observation  at  such  a  point;  therefore  a  new  set  of  data  was  generated  to 
test  the  interpolation.  The  data  from  the  Kriging  was  generated  over  the  range  of  data 
that  was  used  for  generation.  The  results  of  a  comparison  to  NIST  data  can  be  found  in 
Table  15.  The  n  column  is  the  sample  size  and  each  other  number  is  the  percent  error. 


Table  15:  Kriging  Interpolation  Percent  Error  As  Compared  to  NIST  -  Oxygen. 


Data 

n 

H 

Density 

K 

Cv 

cP 

Liquid 

38 

0.28 

0.03 

0.05 

0.07 

0.17 

Sup  Crit 

61 

2.48 

1.95 

0.81 

0.58 

1.55 

Combined 

99 

1.64 

1.21 

0.52 

0.39 

1.02 

Data 

n 

Vise 

U 

S 

Gamma 

Liquid 

38 

1.17 

0.05 

0.13 

0.22 

Sup  Crit 

61 

2.34 

2.18 

0.15 

2.47 

Combined 

99 

1.89 

1.36 

0.14 

1.61 

The  average  errors  for  the  thermodynamic  properties  are  reasonable  when 
averaged  across  the  range  of  states  that  are  of  interest  to  the  rocket  problem.  Error  still 
spikes  near  the  critical  point  as  demonstrated  by  Figure  11. 
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Figure  11:  Kriging  Oxygen  Cp  Error  as  a  function  of  Temperature  and  Pressure. 


Figure  1 1  depicts  a  spike  in  error  of  approximately  20%  in  the  region  of  the 
critical  point  (T  =  278  R,  P  =  73 1  psia).  Based  on  the  method  of  data  collection  for  the 
error  analysis  for  the  super  critical  region,  this  is  a  high  end  estimate  of  error.  The  data 
included  in  the  NPSS  interpolation  tables  for  this  region  is  directly  from  NIST  and 
therefore  any  error  associated  with  it  will  be  from  linear  interpolation  between  two 
points.  For  this  region  in  the  error  analysis,  the  closest  tabulated  value  was  used,  and 
therefore  it  would  be  expected  error  would  be  less  for  a  linear  interpolation.  The 
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approximately  20%  error  in  this  region  is  therefore  a  high  end  estimate  over  a  relatively 
small  region. 

Methane  FPT  Error  Analysis 

Error  associated  with  the  methane  fluid  property  tables  is  significant,  but  was 
controlled  through  the  iterative  improvement  of  the  methods  of  estimation.  As  discussed 
previously,  the  source  NIST  data  is  organized  by  temperature  and  pressure  while  NPSS 
data  must  be  organized  primarily  by  enthalpy  and  pressure.  A  Kriging  estimator 
controlled  through  ModelCenter  reorganized  the  data  via  Kriging  regression  and  created 
the  necessary  tables,  and  NPSS  then  utilized  the  tables  for  property  estimates  through 
linear  interpolation.  The  two  different  interpolation  methods  have  different  methods  of 
reducing  error.  The  Kriging  works  best  when  the  data  is  spread  evenly  over  the  area  of 
interest  with  enough  fidelity  (small  step  size)  to  capture  non-linear  behavior,  but  if  the 
sample  size  is  too  large  (over  1000  points),  the  tool  will  break  down.  On  the  other  hand, 
a  small  step  size  with  data  concentrated  in  the  areas  of  interest  will  reduce  error  for  the 
linear  interpolation. 

The  first  attempt  of  table  generation  was  a  failure  because  the  span  of  data  was 
too  large.  The  problem  was  then  broken  down  by  phase  so  that  two  Kriging  estimators 
could  better  estimate  over  the  range  of  interest.  Upon  examination  of  the  input  data,  it 
was  noticed  the  low  end  of  the  temperature/enthalpy  range  was  neglected  with  regards  to 
data  count.  This  was  remedied  with  the  additional  more  NIST  data.  The  rebuilt  Kriging 
estimators  improved  to  a  useful  point.  The  errors  associated  with  the  Kriging  estimators 
are  in  Table  16. 
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Table  16:  Kriging  Interpolation  Percent  Error  As  Compared  to  NIST  -  Methane. 


n 

H 

Density 

K 

Cv 

Vise 

U 

S 

Gamma 

cP 

153 

9.65% 

1.33% 

1.20% 

0.26% 

2.27% 

4.61% 

11.15% 

2.98% 

3.02% 

Although  the  error  for  enthalpy  and  entropy  is  significant,  in  the  interest  of  time 


and  because  the  error  is  concentrated  in  one  location,  it  was  decided  to  discontinue  the 
effort  of  reducing  the  estimation  error.  This  decision  was  made  because  the  primary 
source  of  estimation  error  remains  at  approximately  20%  in  the  estimation  of  the  heat 
transfer  coefficients  via  the  Bartz  empirical  relationship.  Additionally,  the  error 
consistently  occurs  at  low  temperatures  (200  Rankin  and  below).  The  only  place  this 
state  exists  in  the  model  is  before  the  fuel  pumps.  Because  the  cooling  volume  element, 
CoolingVolume02.int,  uses  enthalpy  and  pressure  to  define  the  state,  the  error  at  this 
location  could  propagate  to  the  other  fluid  properties.  It  is  therefore  recommended  that 
the  low  temperature  methane  error  associated  with  the  Kriging  estimator  be  reduced 
further  before  additional  research  is  undertaken  with  the  methane  fluid  property  table 
files. 

To  demonstrate  that  the  calculation  error  for  enthalpy  is  located  at  the  lower 
temperatures,  the  error  calculation  is  repeated  on  the  same  data  set  after  eliminating  data 
below  200  Rankin.  This  can  be  seen  in  Table  17. 


Table  17:  Kriging  Interpolation  %  Error  Compared  to  NIST,  Methane:  >  200R. 


n 

H 

Density 

K 

cv 

Vise 

U 

S 

Gamma 

cP 

137 

1.03% 

1.44% 

1.23% 

0.24% 

1.96% 

3.12% 

10.27% 

3.26% 

3.36% 

It  is  unreasonable  to  ignore  this  error  because  the  fuel  feed  system  does  have  a 
section  that  is  at  a  low  temperature  and  at  least  one  element  in  this  region  uses  the 
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enthalpy  to  define  the  state  of  the  fluid.  The  large  entropy  error  is  not  a  concern  for  the 
MDEAN  as  none  of  the  MDEAN  elements  employed  entropy  to  define  the  state. 

Combustion  Products  FPT  Error  Analysis 

The  first  several  iterations  of  the  combustion  products  fluid  property  tables  were 
unsatisfactory  with  regards  to  error.  The  error  in  this  table  came  from  two  different 
sources,  the  NPSS  linear  interpolation  and  the  method  of  generation  of  the  reverse  lookup 
table,  temperature  as  a  function  of  enthalpy,  pressure,  and  O/F  ratio,  T=f(H,P,0/F). 

Initial  iterations  had  significant  linear  interpolation  error  due  to  poor  choices  in  step  sizes 
for  the  independent  variables  in  the  regions  of  interest.  Of  the  few  initial  points  of  NPSS 
convergence,  the  area  of  interest  for  this  table  is  defined  as  pressure  from  1600  to  1700 
psia  and  a  temperature  of  6000  to  7000  Rankin.  Nonlinear  behavior  is  detrimental  to  the 
linear  interpolation  process,  but  it  can  be  minimized  through  decreasing  step  sizes  in 
important  regions. 

Initial  guesses  of  step  sizes  in  the  area  of  interest  sometimes  caused  linear 
interpolation  error  on  the  order  of  36%  (at  a  single  point  in  the  area  of  interest).  The 
method  of  reducing  this  error  was  to  add  significantly  more  data  in  the  region.  Increasing 
the  number  of  data  points  from  2500  to  5500  reduced  interpolation  error  considerably. 
The  average  error  across  the  all  of  the  thermodynamic  properties  in  the  area  of  interest  is 
negligible  as  described  in  Table  18. 
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Table  18:  Average  Error  in  Area  of  Interest  by  Thermodynamic  Property. 


Property 

Average  Error 

N 

Enthalpy 

0.00012% 

33 

Density 

0.00005% 

33 

Gamma 

0.00000% 

33 

Cp 

0.000075% 

33 

Thermal  Conductivity 

0.00005% 

33 

Viscosity 

0.00015% 

33 

Entropy 

0.000009% 

33 

Internal  Energy 

0.00046% 

33 

There  is  error  associated  with  gamma,  but  it  is  on  the  order  of  rounding  error  and 
therefore  was  not  seen  in  the  error  calculation.  The  percent  error  was  calculated  by 
comparing  the  NPSS  output  value  of  the  thermodynamic  property  to  the  CEA  calculated 
value.  The  NPSS  output  values  were  obtained  by  exercising  the  solver  over  one  iteration 
and  observing  the  properties  assigned  to  the  exit  port  of  the  RocketCombl  Element.  O/F 
was  held  constant  at  3  for  the  purpose  of  this  measurement.  The  CEA  comparison  data 
was  developed  by  using  the  ModelCenter  carpet  plot  function  on  a  quick-wrapped  CEA. 
Further  error  analysis  would  improve  confidence,  but  because  of  the  small  error 
associated  with  the  linear  interpolation,  further  error  analysis  on  this  table  will  not  be 
undertaken.  Additionally,  the  method  by  which  the  data  was  generated,  through  CEA, 
and  that  the  table  values  are  considered  true  (as  compared  to  the  double  interpolated 
values  of  the  reactant’s  tables)  suggests  that  this  portion  of  the  combustion  products  table 
is  of  comparably  high  quality. 

Another  source  of  error  with  respect  to  the  combustion  products  thennodynamic 
properties  is  the  reverse  lookup  table,  temperature  as  a  function  of  enthalpy,  pressure,  and 
O/F.  This  table  is  notably  employed  in  the  combustion  chamber  to  set  the  state  of  the 
flow  at  the  injector  face.  This  state  and  the  combustion  chamber  exit  port  state  are  used 


83 


to  calculate  the  energy  balance  across  the  combustion  chamber,  which  is  an  important 
dependent  relationship  for  the  MDEAN  convergence.  The  error  associated  with  this  table 
was  originally  due  to  two  interpolations,  first  a  Kriging  to  build  the  table  and  then  a  linear 
interpolation.  The  error  was  significant  and  challenging  to  measure  and  ranged  from  an 
average  of  10%  to  an  average  of  35%  over  the  various  experiments  that  utilized  this 
method.  A  different  approach  was  attempted  that  eliminated  both  sources  of  error  and 
introduced  another  that  was  easily  quantified.  The  fluid  property  tables  can  use  functions 
instead  of  tables  and  this  method  was  used  in  previous  tables  with  only  two  independent 
variables.  The  three  dimensional  function  was  developed  using  a  stepwise,  cubic 
regression  in  ModelCenter.  The  increase  in  data  density  for  the  combustion  products 
tables  discussed  above  provided  enough  data  to  create  a  polynomial  with  suitable 
accuracy  over  the  range  of  possible  enthalpies,  pressures,  and  OF  ratios.  The  adjusted  R~ 
for  the  polynomial  is  96.08%.  Although  this  is  not  a  great  fit  when  compared  to  the 
Kriging  fits  for  the  propellants,  further  analysis  suggests  that  it  is  sufficient  for  this  effort. 

The  error  associated  with  the  polynomial  is  easy  to  calculate  using  a  spreadsheet 
and  the  original  data  used  for  the  polynomial  generation.  The  average  error  associated 
with  the  polynomial  over  the  span  of  data  used  for  generation  is  2.98%,  examined  over 
5500  data  points.  Of  the  530  data  points  that  were  within  the  area  of  interest,  the  average 
error  is  1.22%.  Error!  Reference  source  not  found.There  is  significant  error  when 
using  this  polynomial  model.  However,  most  of  the  high  values  of  error  are  found  at  the 
extremities  of  the  problem,  at  very  low  pressures  or  very  low  enthalpy.  These  conditions 
are  not  usually  associated  with  combustion  chambers.  In  the  region  of  interest,  there  is 
still  error  associated  with  the  polynomial  approximation  of  temperature,  but  the 
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maximum  error  is  about  4%  and  is  associated  with  an  0/F  ratio  of  4.25.  Based  on  earlier 


discussion  of  the  optimal  O/F  ratio  given  an  expansion  ratio,  it  is  unlikely  the  MDEAN 
will  have  optimal  performance  beyond  an  O/F  of  four. 

Another  solution  to  calculate  the  temperature  exists  and  initial  tests  show  that  it  is 
compatible  with  NPSS  fluid  property  table  syntax.  It  can  be  inferred  from  observing 
Error!  Reference  source  not  found,  that  a  better  fit  could  be  obtained  for  the  data  by 
running  cubic  regressions  over  data  sets  of  constant  OF.  Each  OF  ratio  has  a  different 
error  curve,  suggesting  different  polynomials  will  fit  the  different  curve  better  than  a 
unified  polynomial.  Polynomial  development  supports  this  assumption  by  examination 
of  adjusted  R2  values  which  can  be  found  in  Table  19. 


Table  19:  Adjusted  R2  for  Constant  OF,  T  =  f  (H,  P)  Polynomials. 


OF  Ratio 

Adjusted  R2 

#  Data  Points 

2 

96.65 

1250 

3 

95.71 

1250 

3.5 

95.56 

1250 

4 

96.27 

1250 

5 

95.96 

1250 

These  models  obtained  similar  levels  of  successful  fit  with  roughly  20%  of  the 
data  used  in  the  variable  O/F  polynomial.  With  additional  data  and  additional  logic  in  the 
fluid  property  table  function,  an  increase  in  accuracy  might  be  realizable.  The  logic  in 
the  fluid  property  table  is  complex  as  the  tables  do  not  support  simple  if,  then,  else  logic 
operators.  There  is  a  workaround  with  a  C/C++  computational  equivalent,  the  ternary 
operator,  '?’.  In  addition  to  calling  the  proper  functions,  this  method  would  also  require 
an  interpolation  to  use  O/F  ratios  other  than  those  supplied  as  polynomials.  Because  of 
the  added  complexity  of  writing  the  logic  to  include  interpolation,  the  more  straight 
forward  solution,  a  three  dimensional  polynomial,  is  employed. 
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Run  Point  Generation  -  25,000  lbf 

The  model  design  space  was  developed  over  the  five  variables  expansion  ratio, 
throat  area,  chamber  length,  characteristic  length,  and  oxidizer  to  fuel  ratio.  The  first 
exploration  maintained  the  thrust  variable  at  25,000  lbf,  and  it  included  2274  runs  of 
which  541  converged  (23.8%).  The  thrust  level  was  chosen  to  accommodate  the  NGE 
goals  outlined  previously.  The  first  exploration  was  over  the  five  listed  design  variables 
and  included  the  search  outlined  in  Table  20. 


Table  20:  Initial  Search.  25,000  lbf  Thrust. 


Design  Variable 

Low  Value 

High  Value 

Expansion  Ratio 

4 

40 

Throat  Area  (inz) 

8 

20 

Chamber  Length  (in) 

8 

32 

Characteristic  Length  (in) 

40 

61 

O/F  Ratio 

2.2 

4.9 

Number  of  Levels 

4 

Subsequent  runs  expanded  on  areas  of  convergence  and  extended  the  expansion 
ratio  to  60.  The  goal  of  the  high  expansion  ratio  runs  was  to  find  designs  that  met  the 
goal  for  Isp.  In  the  high  expansion  ratio  exploration,  the  rate  of  convergence  dropped  to 
approximately  13%.  The  two  figures  of  merit  calculated  by  the  model  are  specific 
impulse  and  the  thrust  to  weight  ratio.  Figure  12  shows  the  results  of  the  541  converged 
cases  plotted  as  Isp  against  thrust-to-weight.  The  vertical  line  represents  the  NGE  thrust- 
to-weight  goal  of  108  and  the  horizontal  limit  for  Isp,  the  top  of  the  graph,  is  383  seconds 
because  no  points  met  or  exceeded  the  NGE  goal. 
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Figure  12:  Isp  Versus  Thrust-to-W eight.  25,000  lbf  of  Thrust. 

Because  no  points  reach  the  goal  of  383  seconds  of  Isp,  a  point  with  thrust  to 
weight  ratio  greater  than  108  was  initially  selected  for  further  analysis.  The  process  for 
selection  was  based  on  the  concept  of  a  Pareto  frontier,  which  is  described  as  the  subset 
of  values  which  are  Pareto  efficient.  A  Pareto  efficient  point  is  one  where  there  is  no 
opportunity  for  a  Pareto  improvement.  A  Pareto  improvement  is  a  change  from  one  point 
to  another  where  there  is  an  improvement  in  one  value  without  a  corresponding  reduction 
in  the  other  value.  The  area  of  the  above  Figure  12  where  thrust-to-weight  is  greater  than 
108  is  expanded  in  Figure  13.  The  Pareto  efficient  points  are  circled  and  can  be  grouped 
into  two  basic  groups  based  on  Isp  where  the  3  points  on  the  left  are  at  about  350  seconds 
and  the  3  points  on  the  right  are  at  about  340  seconds.  The  point  with  the  highest  thrust- 
to-weight  from  the  grouping  of  the  higher  Isp  values  was  selected  for  further  exploration 
in  the  following  section. 
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Figure  13:  High  Thrust  to  Weight  Ratio  Pareto  Front  at  25,000  lbf  Thrust. 

The  first  point  to  be  analyzed  further  is  described  in  Table  2 1 .  The  initial  point  is 
from  the  data  set  above,  marked  on  Figure  13,  and  the  second  point  is  the  result  of 
parametric  sweeps  over  each  input  design  variable.  The  goal  of  these  parametric  sweeps 
was  a  crude  attempt  at  optimization.  ModelCenter  supports  a  more  thorough 
optimization  package,  but  it  was  not  employed  due  to  time  constraints  and  fragility  of  the 
MDEAN  model.  The  design  point  after  the  parametric  sweeps  serves  as  the  engine 
concept  selected  for  further  developed. 


Table  21:  Design  Point  of  Interest. 


Initial  Point 

After  Parametric  Sweeps 

Expansion  Ratio 

7 

7 

Throat  Area  (in2) 

8.5 

8.4 

Chamber  Length  (in) 

11 

11 

Characteristic  Length  (in) 

50 

50 

O/F  Ratio 

3.3 

3.25 

Specific  Impulse  (s) 

349 

349.3 

Thrust  to  Weight  Ratio 

120 

120.7 
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Alternate  Engine  Concepts 


Because  of  the  automation  capabilities  of  ModelCenter,  the  MDEAN  search 
expanded  into  other  interesting  ranges  of  exploration.  As  discussed  above,  the  first 
exploration  away  from  the  areas  of  likely  convergence  was  high  expansion  ratios.  The 
operational  upper  stage  engines  with  large  values  of  Isp  utilize  bell  nozzles  with  very 
large  expansion  ratios.  For  an  aerospike,  a  large  expansion  ratio  is  directly  related  to  the 
size  of  the  engine.  It  is  expected  that  a  higher  expansion  ratio  (needed  to  achieve  a  higher 
Isp)  is  defined  as  a  larger  engine  and  would  in  turn  have  a  reduced  thrust-to-weight.  Table 
22  describes  one  point  that  reached  the  Isp  goal  after  some  additional  exploration  with 
parametric  sweeps. 


Table  22:  High  Isp  Design  Point  of  Interest  at  25k  lbf  thrust. 


Initial  Point 

After  Parametric  Sweeps 

Expansion  Ratio 

60 

60 

Throat  Area  (in2) 

8 

8 

Chamber  Length  (in) 

14 

14 

Characteristic  Length  (in) 

50 

48 

OF  Ratio 

2.75 

2.75 

Specific  Impulse  (s) 

382.9 

383 

Thrust  to  Weight  Ratio 

13.8 

14.3 

This  point  will  not  be  explored  in  detail  because  of  the  poor  thrust-to-weight  ratio. 
It  is  an  interesting  data  point  because  it  is  representative  of  the  magnitude  of  tradeoff  it 
would  take  to  achieve  the  derived  NGE  Isp  goal.  The  relationship  between  mass, 
expansion  ratio,  and  Isp  is  explored  in  Figure  14.  This  plot  was  generated  by  a  parametric 
sweep  over  expansion  ratio  keeping  all  other  design  variables  constant.  The  static  design 
variables  are  outlined  in  Table  23. 
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Table  23:  Additional  Design  Variables  for  Expansion  Ratio  Sweep  in  Figure  14. 


Throat  Area  (in2) 

9 

Chamber  Length  (in) 

12 

Characteristic  Length  (in) 

50 

Oxidizer  to  Fuel  Ratio 

3 

The  Isp  performance  of  the  rocket  increases  asymptotically  to  a  level  near  the  goal 
of  383  seconds.  At  the  same  time,  the  engine  mass  increases  almost  linearly  so  there  are 
diminishing  returns  to  increasing  the  expansion  ratio  on  the  MDEAN. 
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Figure  14:  Effect  of  Expansion  Ratio  on  Isp,  T/W,  and  Mengine. 
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Alternate  levels  of  thrust  were  explored  in  an  effort  to  find  interesting 
performance  estimates.  First,  thrusts  of  10k,  15k,  and  20k  lbf  were  explored  and  the 
results  are  graphed  as  Isp  versus  thrust-to-weight  in  Figure  15.  None  of  the  explored, 
converged  points  met  the  goal  for  Isp  and  only  one  point  approached  the  goal  for  thrust- 
to-weight.  The  low  thrust  options  were  explored  for  the  possibility  of  combining 
multiple  engines. 
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Figure  15:  Isp  Versus  Thrust-to-Weight.  Thrust  =  10k,  15k,  and  20k  lbf. 


Next,  the  thrust  levels  were  increased  to  30,000  pounds  of  thrust  and  the 
performance  results  can  be  seen  in  Figure  16.  Again,  perfonnance  is  similar  to  that  found 
at  25,000  pounds  of  thrust. 
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Figure  16:  Isp  Versus  Thrust-to-Weight.  Thrust  =  30,000  lbf. 

Figure  17  shows  the  performance  of  the  converged  MDEAN  designs  at  35,000 

pounds  of  thrust.  This  scan  was  over  1296  runs  with  75  design  points  converged  (5.8%). 
This  thrust  level  provides  no  points  that  are  more  interesting  than  the  25,000  lbf  of  thrust 
exploration.  The  lower  rate  of  convergence  suggests  the  cycle  balances  less  often. 
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Figure  17:  Isp  Versus  Thrust-to-Weight.  Thrust  =  35,000  lbf. 
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Exploration  of  Selected  Run  Point 

Several  physical  constraints  are  checked  throughout  the  model  in  an  effort  to 
eliminate  run  points  which  represents  physically  unrealizable  designs.  This  section 
discusses  those  considerations  and  explores  the  feasibility  of  the  chosen  input  parameters. 

The  checks  related  to  the  combustion  chamber  are  related  to  the  gas  flow  velocity 
in  the  combustion  zone  and  the  injection  ratio.  The  Chamber  Mach  number  is  0. 13.  This 
satisfies  the  requirement  that  the  chamber  fluid  velocity  be  “slow.”  The  injection  ratio  is 
the  chamber  flow  area  at  the  injector  divided  by  the  area  at  the  throat  and  should  be 
greater  than  about  2  as  a  rule  of  thumb.  The  Injection  Ratio  is  5.05. 

The  cooling  jackets  are  critical  to  the  convergence  of  the  NPSS  simulation  as  they 
provide  the  energy  flow  for  the  dual  expander  cycle.  The  first  of  the  constraints 
concerning  the  two  cooling  jackets  is  the  fluid  flow  must  not  approach  or  surpass  sonic 
velocity.  The  maximum  fuel  Mach  number  is  0.09283  and  the  maximum  oxidizer  Mach 
number  is  0. 15837.  These  seem  acceptable  if  not  a  bit  slow.  The  coolant  velocity  in  the 
cooling  channels  is  related  to  the  rate  of  heat  flow  across  the  wall.  This  suggests  there  is 
plenty  of  margin  to  reduce  the  depth  of  the  channels  to  increase  the  velocity  of  the  fluids. 
If  the  wall  temperatures  of  the  cooling  jackets  are  too  hot,  it  would  be  important  to 
increase  the  fluid  velocity.  The  wall  temperature  profile  is  shown  in  Figure  18.  The 
breakdown  point  of  silicon  carbide  is  5405  Rankin  and  it  is  shown  as  the  horizontal  line 
on  Figure  18.  The  X  axis  on  Figure  18  is  axial  distance  from  the  injector  face  in  inches 
with  the  throat  at  1 1  inches.  The  two  final  stations  on  the  spike  side  cooling  jacket  are  on 
the  aerospike,  outside  of  the  engine.  The  maximum  temperature  is  on  the  spike  side  at  a 
temperature  of  353 1  Rankin.  This  is  a  significantly  higher  wall  temperature  than  was 
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experienced  in  Hall’s  DEAN  model,  about  1000  Rankin  [4].  One  advantage  of  the 
cryogenic  hydrogen  is  its  excellent  coolant  qualities.  This  difference  in  cooling 
effectiveness  necessitated  a  change  in  the  spike  cooling  jacket  material  from  oxygen  free 
copper  to  silicon  carbide.  Silicon  carbide  has  favorable  properties  at  high  temperatures 


but  creates  challenges  in  manufacturing. 


Axial  Distance  from  Injector  (in) 


Spike  Side  Wall 
Temperature 


SiC  Decomposition 
(5405R) 


Chamber  Side  Wall 
Temperature 


Figure  18:  Wall  Temperature  Profiles  for  Cooling  Jackets. 


In  addition  to  the  temperature  profile,  the  pressure  profile  is  important  to  the 
proper  function  of  the  cooling  cycle.  Figure  19  shows  the  pressure  drops  across  the  two 
cooling  systems.  The  curves  seem  reasonable  and  the  fuel  pumps  increase  the  fuel 
pressure  by  equal  amounts,  as  expected  in  this  version  of  the  MDEAN. 
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Figure  19:  Propellant  Feed  System  Pressure  Profile. 


An  important  constraint  on  the  pressure  profile  is  the  injector  pressure  drop  which 
should  be  approximately  20%  to  50%  of  the  chamber  pressure.  The  oxidizer  pressure 
drop  across  the  injector  is  3 17  psia,  or  19.56%  of  chamber  pressure.  The  fuel  pressure 
drop  across  the  injector  is  375  psia,  or  23%  of  chamber  pressure.  Although  the  oxidizer 
side  is  a  bit  low,  this  was  deemed  acceptable.  It  would  be  ideal  to  increase  the  pressure 
of  the  fluid  before  the  injector.  As  it  stands,  there  is  insufficient  margin  to  increase  the 
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fluid  velocity  through  the  oxidizer  cooling  jacket.  Although  the  wall  temperature  could 
be  reduced,  the  pressure  drop  through  the  cooling  jacket  would  increase,  reducing  the 
pressure  at  the  injector. 

The  final  check  is  against  the  aerospike  thickness.  The  walls  of  the  aerospike  are 
built  up  by  three  layers.  The  cooling  jacket  is  the  interface  with  the  combustion  side  hot 
gases.  Below  the  cooling  jacket  are  the  cooling  channels  of  defined  thickness  and  depth. 
Supporting  these  is  a  structural  jacket.  The  remainder  of  the  radius  (near  the  center  of  the 
spike)  is  hollow.  The  thickness  of  each  component  and  the  radius  of  the  spike  at  every 
station  are  all  independently  calculated.  For  every  station,  the  three  layers  of  the  spike 
are  checked  against  the  radius  of  the  spike.  If  the  sum  of  the  layers  is  greater  than  the 
spike,  the  MDEAN  rejects  the  run  point. 

The  physical  dimensions  of  the  MDEAN  are  displayed  in  Figure  20.  The  small 
physical  size  of  the  DEAN  family  of  rockets  is  an  important  characteristic. 


Figure  20:  The  MDEAN  Physical  Dimensions. 
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The  MDEAN  volume  is  well  under  the  NGE  spatial  requirements  of  90  inches 
from  gimbal  to  nozzle  exit  and  an  exit  diameter  of  73  inches  [5].  For  those  missions  that 
have  challenging  volume  requirements,  the  MDEAN ’s  comparatively  small  size  might  be 
beneficial.  A  scale  comparison  of  the  NGE  requirements  to  the  selected  MDEAN  engine 
is  shown  in  Figure  21.  Additional  length  savings  can  be  realized  through  truncation  of 
the  aerospike  with  a  small  loss  of  perfonnance  [8,  13]. 


Figure  21:  MDEAN  Overlaying  the  NGE  Spatial  Requirements. 

Comparison  of  MDEAN 

In  this  section,  the  MDEAN  is  compared  to  several  other  rocket  concepts.  Table 
24  shows  the  primary  perfonnance  parameters  for  several  different  rockets,  now 
including  the  MDEAN. 
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Table  24:  Comparable  Engine  Performance. 


Engine 

MDEAN 

RD-185 

RD-167 

RD-160 

NGEDer 

NGE 

H-DEAN 

Fuel 

ch4 

ch4 

ch4 

ch4 

ch4 

h2 

h2 

Thrust  (lbf) 

25K 

40,344 

79,366 

4,409 

25K 

25K 

44,694 

Up  (s) 

349.3 

378 

379 

380.6 

383 

465 

429.8 

T/W 

120.7 

44 

63 

15.5 

108 

N/A 

142.2 

Mass  (lbm) 

207 

917 

1260 

284 

231 

N/A 

314 

Mprop  (lbm) 

22,323 

17,715 

17,588 

17,388 

N/A 

N/A 

20,600 

The  mass  of  propellant  for  each  rocket  was  calculated  for  the  Av  calculated  in 
chapter  2,  4.3  km/s.  The  payload  mass  is  assumed  as  5000  kg.  The  finert  estimated  for  the 
MDEAN  was  extended  to  the  other  methane  rockets.  As  expected,  the  MDEAN  requires 
the  greatest  amount  of  propellant  as  it  has  the  lowest  Isp.  For  the  MDEAN  to  have  an 
advantage  over  the  other  rockets,  the  high  thrust-to-weight  ratio  must  become  an 
advantage.  The  weight  savings  on  hardware  must  be  similar  to  the  weight  losses  due  to 
inferior  propellant  efficiency.  Figure  22  shows  the  mass  of  required  propellant  versus  the 
change  in  velocity  for  the  methane  engines. 
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Figure  22:  Methane  Rockets  Propellant  Mass  Vs.  Av. 
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It  is  clear  that  for  any  maneuver  of  significant  Av,  Isp  is  the  more  important  figure 
of  merit.  The  MDEAN  should  have  a  perfonnance  advantage  over  the  other  methane 
rockets  for  a  required  change  in  velocity  of  less  than  2  km/s. 

One  possible  benefit  of  the  methane  rocket  over  a  similar  hydrogen  rocket  is 
savings  in  the  mass  of  the  propellant  tanks.  This  comparison  is  made  between  Hall’s 
hydrogen  DEAN  [8]  and  the  MDEAN.  The  method  to  estimate  the  mass  of  the 
propellant  tanks  begins  with  the  change  in  velocity  calculation  from  chapter  2.  The 
transfer  orbit  used  for  the  calculation  has  a  Av  requirement  of  4.3  km/s.  The  payload  for 
delivery  to  orbit  is  assumed  to  be  5,000  kg.  Next,  a  finert  is  assumed  for  the  rocket  stage 
that  is  defined  by  equation  52  [4]. 


ft 


m 


inert 


inert 


m-prop  4”  t7l 


inert 


(52) 


where 


minert  =  Vehicle  mass  excluding  propellant  and  payload  (kg) 
mProp  =  Mass  of  required  propellant  (kg) 

The  mass  of  propellant  can  then  be  calculated  by  combining  the  definition  of  finert 
with  the  ideal  rocket  equation  to  reach  equation  53. 


f— ) 

g  VspdoJ 


1(1  -finert) 
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mprop  /  Av  \ 

1  _  f.  pVspdo) 

J  inert c 

With  the  O/F  Ratio,  the  mass  of  each  propellant  can  then  be  calculated.  With  the 
total  mass  of  each  propellant  and  the  desired  thermodynamic  state  of  each  storage  tank, 
the  density  can  be  detennined  through  the  NIST  [25,  26]  and  the  volume  of  the  tanks  can 
be  calculated.  At  this  point,  there  are  several  methods  for  the  calculation  for  the  mass  of 
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cryogenic  fluid  storage  tanks.  The  tank  mass  calculation  was  chosen  based  on  the  rule  of 
thumb  that  the  tank  weight  should  be  about  5-15%  of  the  propellant  weight  [29].  The 
method  of  tank  mass  calculation  is  taken  from  Humble  et  al.  [4],  Once  the  tanks  are 
sized,  a  pressurization  system  mass  is  estimated,  providing  the  mass  of  all  inert 
components.  This  allows  for  the  calculation  of  finert  and  an  iterative  process  to  solve  for 
the  tank  masses.  The  results  of  this  analysis  can  be  found  in  Table  25. 


Table  25:  MDEAN  Compared  to  DEAN. 


MDEAN 

DEAN 

Fuel 

Methane 

Hydrogen 

ITlengine  (kg) 

94 

160 

mDr0D  (kg) 

20,245 

23,078 

mfuei  (kg) 

4,764 

3,297 

m<i\idizer  (kg) 

15,482 

19,781 

Illtanks  (kg) 

792 

2103 

Hlinert  (kg) 

3067 

8,050 

m  pay  load  (kg) 

5,000 

5,000 

mtotai,  (including  payload)  (kg) 

28,522 

35,952 

The  tank  masses  were  based  on  spherical  tanks  of  the  materials  chosen  by  Hall  for 
the  hydrogen  (INCOLOY  909)  and  oxygen  (INCONEL  718)  feed  systems  [8].  The  mass 
was  calculated  by  detennining  the  thickness  of  the  tank  assuming  a  safety  factor  of  2,  or  a 
burst  pressure  of  400  psia.  When  compared  to  empirical  relationships  for  cryogenic 
tanks,  all  of  the  tanks  are  light  except  for  the  hydrogen  tank,  which  is  significantly 
heavier.  This  is  not  a  perfect  comparison  as  the  DEAN  operates  at  50,000  lbf  of  thrust 
while  the  MDEAN  is  designed  for  25,000  lbf. 
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V.  Conclusions  and  Recommendations 


This  chapter  first  explores  the  conclusions  and  significance  of  the  MDEAN 
research.  Then,  based  on  experience  with  the  model,  recommendations  for  future  work 
are  made. 

Conclusions  of  Research 

Based  on  the  runs  presented  in  chapter  4,  the  MDEAN  does  not  represent  a 
significant  boost  over  the  state  of  the  art.  It  does  not  meet  both  requirements  of  383 
seconds  of  specific  impulse  with  a  thrust  to  weight  ratio  of  over  108.  The  estimates  for 
all  of  the  simulated  methane  engines  could  prove  useful  to  an  actual  multi  stage  design. 
The  MDEAN  provides  an  advantage  based  on  volume  savings  when  compared  to  a 
hydrogen  rocket  or  a  rocket  with  a  bell  nozzle.  There  could  be  a  mission  that  would 
favor  the  characteristics  of  a  MDEAN  rocket. 

Significance  of  Research 

The  MDEAN  adds  significant  capability  over  the  existing  DEAN  capability.  The 
fluid  property  tables  can  be  used  in  other  models  with  some  attention  to  the  low 
temperature  methane  error.  Additionally,  this  paper  outlined  the  process  for  exploring 
different  propellants  with  NPSS.  By  referencing  the  methods  used  for  the  development 
of  the  fluid  property  tables,  it  should  be  possible  to  develop  tables  for  alternate  chemistry 
in  a  timely  manner. 

Additionally,  the  exploration  of  an  upper  stage  dual  expander  aerospike  nozzle 
rocket  engine  provided  a  large  set  of  engine  mass/performance  design  data.  Given 
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mission  and  geometric  constraints,  it  should  be  possible  to  query  this  research  to 
detennine  if  there  is  a  liquid  methane/liquid  oxygen  design  concept  that  would  provide  a 
reasonable  solution.  At  the  least,  it  should  be  possible  to  deconstruct  requirements  into  a 
range  for  the  design  input  variables  and  explore  the  solution  space. 

Recommendations  for  Future  Research 

The  first  several  recommendations  involve  improvements  to  the  existing  model. 
First,  the  methane  Kriging  estimator  for  the  fluid  property  tables  requires  an 
improvement.  The  error  at  low  temperatures  should  be  reduced  by  including  more  low 
temperature  data  and  rebuilding  the  estimator  in  ModelCenter.  This  is  a  relatively 
straight  forward  process,  but  was  not  accomplished  due  to  time  constraints.  Next,  the 
MDEAN  should  be  run  over  a  large  set  of  design  input.  Because  of  time  constraints,  the 
model  was  mostly  explored  over  the  six  input  design  variables:  expansion  ratio,  throat 
area,  chamber  length,  characteristic  length,  oxidizer-to-fuel  ratio,  and  thrust.  This  design 
space  could  be  opened  up  to  look  at  the  aspect  ratio  of  the  cooling  channels,  the  number 
of  cooling  channels,  the  percent  of  oxidizer  flow  through  the  bypass  and  more. 

Before  the  design  space  is  fully  explored,  an  upgrade  to  decrease  the  time  per  run 
of  the  model  should  be  explored.  One  opportunity  to  decrease  the  time  per  run  is  to  leave 
NPSS  open  from  run  to  run.  Currently,  for  each  ModelCenter  run,  NPSS  is  opened,  ran, 
and  then  closed.  By  just  leaving  NPSS  open  from  run  to  run  could  save  significant  time 
over  thousands  of  runs.  Averaging  over  50  runs,  the  model  takes  approximately  8 
seconds  to  complete.  Some  failure  modes  within  NPSS  are  on  the  order  of  two  minutes 
(not  seen  in  the  measured  50  runs),  but  that  is  usually  because  the  design  point  failed  to 
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converge  within  the  limit  of  iterations.  This  could  and  should  be  improved  by  reducing 
the  maximum  number  of  iterations  from  5000. 

A  test  plan  for  the  model  was  developed,  but  not  fully  executed.  More  time 
should  be  spent  on  thorough  validation  and  verification.  Much  of  the  work  on  the 
MDEAN  and  the  concurrently  developed  DEAN7  was  spent  on  debugging  reused  code. 

More  data  could  be  used  in  several  parts  of  the  model.  More  high  temperature 
methane  and  oxygen  data  would  provide  a  better  buffer  for  the  cooling  channels. 
Currently,  the  state  of  the  fluid  in  the  cooling  channels  can  be  at  a  higher  temperature 
than  is  included  in  the  fluid  property  tables  or  the  sonic  velocity  constraint  checkers. 

This  is  because  of  a  limitation  in  the  available  NIST  data.  More  data  would  prevent  this 
model  failure  mode,  and  if  a  source  of  thermochemical  data  is  identified,  more  propellant 
fluid  property  tables  could  be  created.  Additionally,  more  material  property  data  as  a 
function  of  temperature  is  needed  for  silicon  carbide.  The  model  is  currently  built  around 
two  data  points  for  silicon  carbide  ultimate  tensile  strength  and  yield  strength,  at  540  R 
and  1080  R. 

NPSS  has  been  advanced  to  a  commercial  version  of  2.3  while  DEAN7  and 
MDEAN  are  built  with  NPSS  version  1.65.  All  NPSS  rocket  elements  would  need  an 
upgrade  into  2.3,  but  they  should  somewhat  similar  because  both  versions  use  a  syntax 
derived  from  C++.  Additionally,  the  rocket  elements  could  possibly  be  upgraded  to 
include  reacted  fluid  stations  (or  new  rocket  elements  could  be  developed).  DEAN  only 
uses  unreacted  fluid  stations  and  therefore  pushed  the  MDEAN  thermochemical  solution 
into  using  fluid  property  tables.  Several  other  solution  methods  exist  with  different  fluid 
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station  types.  For  example,  NPSS  can  call  CEA  directly  to  calculate  fluid  properties  of 
combustion  products,  but  not  with  unreacted  fluid  stations. 

Summary 

Based  on  several  iterations  of  a  liquid  hydrogen  and  liquid  oxygen  model  built  up 
over  years  of  work,  the  MDEAN  furthers  the  capability  to  analyze  a  diverse  selection  of 
rocket  concepts.  While  failing  to  meet  all  perfonnance  goals  derived  from  the  NGE 
solicitation,  the  MDEAN  yet  provides  relevant  perfonnance  estimates  to  interesting 
rocket  concepts.  The  final  rocket  engine  concept  has  an  Isp  of  349.3  seconds,  a  thrust-to- 
weight  ratio  of  120.7,  and  a  total  engine  mass  of  207.2  lbm.  An  examination  of  this  run 
point  demonstrated  it  is  physically  realizable. 
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Appendix  A:  Simple.mdl 


Due  to  the  complexity  of  the  DEAN  model,  it  was  undesirable  to  experiment  with 
changing  the  fluid  property  tables.  In  addition  to  updating  the  specified  thermodynamics, 
the  fuel  pump  and  the  initial  guesses  of  the  cooling  channel  fluid  states  required  update. 
To  reduce  the  level  of  complexity  of  the  software  update,  a  simple  pressure  fed  rocket 
model  was  developed,  simple.mdl.  Whenever  the  fluid  property  tables  were  updated, 
they  were  first  tested  with  the  simple.mdl  model.  This  allowed  for  an  iterative  software 
update  where  individual  steps  of  the  update  could  be  verified.  The  model  includes 
propellant  tanks,  simplified  frictional  line  losses,  a  rocket  combustion  chamber,  and  a 
rocket  nozzle.  In  addition  to  the  standard  independent  and  dependent  variables  included 
with  the  elements,  the  following  were  added: 


Table  26:  Simple.mdl  Independent  Variables  and  Dependent  Conditions. 


Independent  Variables 

Dependent  Variables 

Throat  Area 

Mass  Flow  Rate 

Chamber  Area 

Thrust 

The  solver  perturbs  the  throat  area  and  the  chamber  area  while  trying  to  attain  the 
specified  mass  flow  rate  and  thrust.  The  development  of  this  model  followed  the  rocket 
design  process  found  in  Space  Propulsion  Analysis  and  Design  [4],  which  provided 
approximate  values  of  all  input  variables  for  the  different  elements.  The  perfonnance  of 
the  model  was  not  examined  as  it  was  simply  a  tool  to  test  the  new  fluid  property  tables 
independently  of  other  model  updates.  The  code  is  presented  below. 


// 

//  This  is  a  simple  rocket  engine. 

//  It  will  serve  as  a  starting  point  for  more  detailed  models. 
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// 

//  01  Nov  11 

// 


MODELNAME  =  "Bipropellant"; 

string  Model  lD  =  "Methane  Oxygen  chemical  engine"; 

0^  Thermochemistry  source 

setThennoPackage("FPT"); 


^/^^c${c${c$!c${c${c${c$!c${c${c${c$}c${c$}c${c$!c${c${c${c$!c${c${c${c$!t${c${c${c  IlTpUt  V  anables 


real  mdot_Ox  =  16.4134;  //  02  fluid  mass  flow  rate  (lbm/s) 

real  mdot_Fu  =  7. 13636;  //  Fuel  fluid  mass  flow  rate  (lbm/s)  ! ! ! ! ! Set  to 

RP- 1  needs  adjustment  in  this  model  to  H2 !!!!!!!!!!!!!! 

real  mdot  =  mdot_Ox  +  mdot_Fu;  //  Total  mass  flow  rate  (lbm/s) 

real  thrust  =  8000.0;  //  target  thrust  (lbf)  7937 

real  P_chamber  =  101.526;  //  chamber  pressure  (lbf/in2  aka  psia) 

//***************************  *  Setup  Elements*  **************************** 


Element  Starter  TankOx  { 


OFR 

=  1.0; 

Pt 

=  140.0; 

//  Pressure  of  ox  tank  given  in  psia 

Tt 

=  162.0; 

//  BP  of  oxygen,  storage  as  liquid  at  high  pressure 

deg  R 

comp 

} 

=  "02NIST"; 

Element  Pipe  OxPipe  { 

Cf 

} 

=  100.0; 

Element  Starter  TankFu  { 

OFR 

=  0; 

Pt 

=  140.0; 

//  Pressure  of  ox  tank  given  in  psia 

Tt 

=  210; 

//  BP  of  CH4  (@latm),  but  storage  as  liquid  at  high 

pressure  deg  R 

comp 

=  "METHANE"; 

//  Tt 

=  36; 

//  BP  of  CH4  (@latm),  but  storage  as  liquid  at  high 

pressure  deg  R 

//  comp  =  "HYDROGEN"; 

} 
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Element  Pipe  FuPipe  { 
Cf  =  100.0; 

} 


Element  RocketCombl  Chamber  { 

OFR  =0.75;  // Oxidizer  fraction  OFR=OF/(l+OF) 

comp  =  "CCH4  02";  //  Composition  of  hot  side  gases,  only  for 

combusting  H2  and  02 

//  comp  =  "C02H2";  //  Composition  of  hot  side  gases,  only  for 

combusting  H2  and  02 

radius_tc  =6.367;  // Radius  of  the  thrust  chamber  just  before 

converging  to  nozzle  (in) 

volume  =1491.5;  //  Volume  of  gas  in  chamber,  in3 

} 


Element  RocketNozzle  Nozz  { 

AR  =100.0;  //  Area  ratio  for  nozzle  exit 

Ath  =  42.09354;  //  Area  of  throat,  given  in  in*  *2 

Ps  =  0.00 1 ;  //  Ambient  pressure  for  upper  stage,  assumed  psia 

realLossCoef  =  1.0;  //  Loss  coefficient  accounting  for  everything  but 

expansion  losses  (set  to  default) 

s_Q  =  1 .0;  //  Scalar  on  heat  transfer  rate  (set  to  default) 


hnkPorts("TankOx.Fl_0",  "OxPipe.Fl  I",  "Ox  piping");  //links  oxidizer 
tank  output  to  a  pipe 

linkPorts("OxPipe.Fl_0",  "Chamber.Floxid",  "Floxidchmbrinput");  //links 
oxidizer  pipe  output  to  chamber  oxidizer  input 

linkPorts("TankFu.Fl_0",  "FuPipe. FI  I",  "Fu  piping");  //links  fuel  tank 
output  to  a  pple 

linkPorts("FuPipe.Fl_0",  "Chamber.Fu  I",  "Fl  fuel  chmbrinput");  //links 
fuel  pipe  output  to  chamber  fuel  input 

linkPorts("Chamber.Fl_tc",  "Nozz.Fl  I",  "Fl  chmbr  to  nozz");  //links 
combustion  chamber  to  the  nozzle 


//* 


Setup  Solver* 


//  Defining  output  variables 
Chamber.Pttc  =  Pchamber; 
Chamber. Fltc.W  =  mdot; 


//Initial  guesses  required  for  solver  initiation 
Chamber.OFR  =  0.75;  //OFR  =  OF  /  (1+OF) 
Chamber.Tt_tc  =  5400;  //(R) 
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TankOx.Pt  =  140.0; 

TankOx.Tt  =  162.0;  //(R) 

TankFu.Pt  =  140.0; 

TankFu.Tt  =210.0;  //(R) 

Chamber.  Winj  =  mdot; 

Chamber. Ful.htRef  =  2391.6892;  //This  is  a  correction  factor  to  get  from  NIST 
to  CEA  Enthalpy  reference  state 

Chamber.Floxid.htRef  =  1 17.0293;  //This  is  a  correction  factor  to  get  from  NIST 
to  CEA  Enthalpy  reference  state 

//Independent  Variables 

Independent  AThroat  { 

varName  =  "Nozz.Ath"; 
auto  Setup  =  TRUE; 


Independent  AChamber  { 

varName  =  "Chamber.radius_tc"; 
auto  Setup  =  TRUE; 


//Dependent  Variables 

Dependent  MassFlowj  //massflow  is  given  in  the  SPAD  description 

eq_lhs  =  "Nozz.W";  //solved  for  in  the  RocketNozzle  Nozz  element 
eq_rhs  =  "Chamber.Winj";  //provided  in  the  problem  statement 
auto  Setup  =  TRUE; 


Dependent  Thrust  { 

eq_lhs  =  "Nozz.Fg"; 
eq_rhs  =  "thrust"; 
auto  Setup  =  TRUE; 


solver.  solutionMode  =  "STEAD  YSTATE"; 
presolverSequence  =  {}; 

autoSolverSetup(); 
solver. maxlterations  =100; 
solver.maxJacobians  =  10; 
run  (); 

cout «  endl; 
cout «  endl; 
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cout «  "  Converged  "  «  solver,  converged; 
cout «  endl; 
cout «  endl; 

cout «  "  ITT  PASS  JAC  BROY  "; 
cout «  endl; 
cout.  width  =  6; 

cout «  solver.iterationCounter 
«  solver.passCounter 
«  solver. numJacobians 
«  solver.numBroydens; 
cout «  endl; 
cout «  endl; 

cout.width  =  4; 

//cout.precision  =  4; 
cout «  "  Oxygen  Tank  "  ; 
cout«"  Pt  "  «  TankOx.Pt; 

cout«"  Tt  "  «  TankOx.Tt; 

cout «  endl; 

cout  «  "  Fuel  Tank  "  ; 
cout  «  "  Pt  "  «  TankFu.Pt; 

cout «  "  Tt  "  «  TankFu.Tt; 

cout «  endl; 

cout.precision  =  4; 

cout.width  =  4; 

cout «  "  Ox  Pipe  " ; 

cout  «  "Pout  "  «  0xPipe.Fl_0.Pt; 

cout «  "  hout  "  «  0xPipe.Fl_0.ht; 

cout «  "  Tout  "  «  0xPipe.Fl_0.Tt; 

cout «  "  rho  out  "  «  OxPipe.Fl_0.rhot; 

cout «  "  comp  out  "  «  OxPipe.Fl_0.comp; 

cout «  "  dp  "  «  OxPipe.dP; 

cout «  "  "; 

cout.width  =  4; 

//cout.precision  =  4; 

cout «  "  W  "  «  OxPipe.Fl  O.W; 

cout «  endl; 

cout «  "  Fu  Pipe  " ; 

cout.width  =  4; 

cout.precision  =  4; 

cout  «  "Pout  "  «  FuPipe.Fl_0.Pt; 

cout «  "  hout  "  «  FuPipe.Fl_0.ht; 

cout «  "  Tout  "  «  FuPipe.Fl_0.Tt; 
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cout «  "  rhoout  "  «  FuPipe.Fl_0.rhot; 

cout «  "  comp  out  "  «  FuPipe.Fl_0.comp; 

cout «  "  dp  "  «  FuPipe.dP; 

cout «  "  "; 

cout.  width  =  3; 

cout.precision  =  3; 

cout «  "  W  "  «  FuPipe.Fl  O.W; 
cout «  endl; 

cout «  "  Chamber  "  ; 
cout.  width  =  4; 
cout.precision  =  4; 

cout «  "  Pt "  «  Chamber.Pttc; 

cout «  "  Tt  "  «  Chamber.Tttc; 

cout.  width  =  4; 
cout.precision  =  4; 

cout«"  Winj  "«  Chamber. Winj; 

cout «  endl; 
cout «  "  "; 

cout «  "  Rtc  "  «  Chamber.radiustc; 

cout «  "  " ; 

cout«"  Wex  "  «  Chamber.Wnozzle; 
cout «  endl; 
cout «  endl; 

cout.  width  =  3; 

cout.precision  =  3; 

cout «  "  Nozzle  " ; 

cout «  "  At  "  «  Nozz.Ath; 

cout.  width  =  5; 

cout.precision  =  5; 

cout «  "  Fg  "  «  Nozz.Fg; 

cout.  width  =  4; 

cout.precision  =  4; 

cout«"  Wnox  "  «  Nozz.Fl  I.W; 
cout «  "  Isp  "  «  Nozz.Isp; 
cout «  endl; 
cout «  endl; 


111 


Appendix  B:  Software  Changes  from  DEAN7  to  MDEAN 

This  appendix  lists  the  changes  that  were  required  to  update  the  standard  DEAN 
model  (version  7)  to  accommodate  methane  fuel  via  fluid  property  tables.  The  first  table 
outlines  the  changes  to  the  NPSS  model.  For  these  changes  to  function  properly,  the 
fluid  property  tables  must  be  placed  in  the  following  directory  in  the  NPSS  folder  which 
is  usually  installed  in  the  program  files  directory  for  a  Windows  install: 

C:\Program  Files  (x86)\NPSS.nt.V165-OPT-Full\DLMComponents\nt 

This  will  of  course  change  based  on  the  chosen  operating  system,  NPSS  version, 
and  file  structure.  This  installation  was  on  Windows  7  with  NPSS  version  1.65.  Note 
that  everything  above  line  68  of  the  NPSS  model  is  also  controlled  through  ModelCenter 
input  variables.  Variables  that  are  considered  design  variables  are  not  included  in  these 
lists.  Changes  are  also  required  at  the  ModelCenter  level  and  are  outlined  in  Table  28. 


Table  27:  NPSS  Updates 


Line 

Update 

Action 

22 

Targetlsp  =  350; 

Update  to  realistic  CH4/02  guessed  Isp 

25 

real  T  TankH  =  170.0; 

Update  to  realistic  CH4  storage  temperature 

39 

real  Nmech  RPMF  =  50000; 

Update  to  fuel  pump  rpm  based  on  maintaining 
specific  speed 

78 

real  c  star  =  6070; 

Update  c  star  value  to  6070  for  methane 

87 

real  ht  OCV1  = -62.312; 

Update  Oxidizer  initial  enthalpy  guesses  to 

NIST  thermodynamic  basis. 

88 

real  ht  OCV2  = -54.0; 

89 

real  ht  CV06  = -52.0; 

90 

real  ht  CVO5  =  -19.0; 

91 

real  ht_CV04  =  5.0; 

92 

real  ht  CV03  =  27.0; 

93 

real  ht  CV02  =  48.0; 

94 

real  ht  CVOl  =  68.0; 

95 

real  ht  OCV3  =  68.0; 

96 

real  ht  TBCV1  =68.0; 

97 

real  ht  OCV4  =  59.0; 

99 

real  ht  HCV1  = -25.243; 

Update  Fuel  initial  enthalpy  guesses  to  methane 
values  from  the  NIST  data. 

100 

real  ht  HCV2  =  -15.6; 
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101 

real  ht  HCV3  =  65.95; 

102 

real  ht  CVH8  =  111.7; 

103 

real  ht  CVH7  =  221.2; 

104 

real  ht  CVH6  =  226.7; 

105 

real  ht  CVH5  =  328.1; 

106 

real  ht  CVH4  =  388.6; 

107 

real  ht  CVH3  =  439.6; 

108 

real  ht  CVH2  =  485.99; 

109 

real  ht  CVH1  =  532.1; 

110 

real  ht  HCV4  =  532.1; 

111 

real  ht  HCV5  =  440.6; 

125 

real  rho  HD1  =  27.89; 

Update  to  initial  methane  guesses  for  density. 

126 

real  rho  DuctH8  =  22.5; 

127 

real  rho  DuctH7  =  19.82; 

128 

real  rho  DuctH6  =  16.116; 

129 

real  rho_DuctH5  =  13.6688; 

130 

real  rho  DuctH4  =  10.45; 

131 

real  rho  DuctH3  =  8.829; 

132 

real  rho  DuctH2  =  7.524; 

133 

real  rho  DuctHl  =  6.649; 

134 

real  rho  TBH0  =  6.66; 

135 

real  rho  HV  =  4.43299; 

537 

setThermoPackage("FPT"); 

Specifies  new  thermodynamic  package 

542 

setDefaultComposition("C  CH4 
02"); 

For  the  chamber  and  nozzle,  the  composition  is 
the  combustion  products  unless  otherwise 
specified.  This  calls  the  fluid  property  table 
named:  C  CH4  02.fpt 

544 

comp  =  "CCH4  02"; 

Specifies  combustion  products  for  combustion 
chamber  exit  port,  redundant  to  line  542,  but 
original  NPSS  model  included  it. 

547 

Ful.comp  =  "METHANE"; 

Specifies  METHANE. fpt  fluid  property  table. 

548 

FI  oxid.comp  =  "02  NIST"; 

Specifies  02_NIST.fpt  fluid  property  table. 

636 

setDefaultComposition("02  NIS 
T"); 

Specifies  02_NIST.fpt  fluid  property  table  for 
all  flow  stations  until  new  composition  is 
specified. 

859 

setDefaultComposition("METHA 

NE"); 

Specifies  ME  THANE,  fpt  fluid  property  table 
for  all  flow  stations  until  new  composition  is 
specified. 

1287 

COMB.Ful.htRef  =  2391.6892; 

Added  correction  factor  to  convert  from  NIST 
thermodynamic  enthalpy  basis  to  CEA  enthalpy 
basis  for  methane. 

1288 

COMB.F1  oxid.htRef  = 

117.0293; 

Added  correction  factor  to  convert  from  NIST 
thennodynamic  enthalpy  basis  to  CEA  enthalpy 
basis  for  oxygen. 
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Table  28:  ModelCenter  Updates. 


Component 

Update 

Action 

DesignV  ariables  .Materi 
als 

Spike  Cool  Mat  =  Silicon 
Carbide 

Change  the  aerospike  side 
cooling  jacket  to  silicon 
carbide  as  the  copper  melts. 

Targetlsp 

Change  to  350 

Redundant  to  NPSS 
change,  but  necessary  if 
NPSS  model  is  to  be  used 
independent  of 

ModelCenter 

DEAN 

TTankHNmechRPMF 

Change  thse  values  to 
match  the  NPSS  input 
values  for  consistency. 

Constraints .  LH2 .  LH2 
Machs. SonVel  H2 

Methane  Kriging  Sonic  Velocity 
estimator 

Replace  H2  SonVel 
estimator  with  CH4  SonVel 
estimator. 

Perfonnance.CEA 

Different  Component  that  selects 
correct  fuel 

Updated  with  different 
wrapper. 
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